Import, scaling & normalisation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing data, scaling and normalisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is intended as an introduction to some of the basic functionality of Jupyter Notebooks used throughout this set of tutorials. First, we demonstrate how to import external packages, read the raw study data and configure the plotting interface. \n",
    "\n",
    "We also demonstrate how to perform Probabilistic Quotient Normalisation, apply different types of scaling to the data matrix and visualise their impact on the raw data matrix and interpretation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Package/code import and environment set up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell imports the various necessary python packages. This cell is run at the beginning of each notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn import preprocessing\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PCA and PLS/PLS-DA/OrthogonalPLS codes provided in these tutorial notebooks are bundled as a python module (package) named **pyChemometrics**. As an example, the next cell imports the *pyChemometrics* scaler and PCA objects. The syntax and functions of these objects will be explained throughout the tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example import\n",
    "from pyChemometrics.ChemometricsPCA import ChemometricsPCA\n",
    "from pyChemometrics.ChemometricsScaler import ChemometricsScaler\n",
    "\n",
    "# Other available objects\n",
    "# 1) ChemometricsPLS\n",
    "# 2) ChemometricsOrthogonalPLS\n",
    "# 3) ChemometricsPLSDA\n",
    "# 4) ChemometricsOrthogonalPLSDA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell sets up the figure display mode. The *notebook* mode allows interactive plotting. Another option is to select *inline*, to obtain static plots in a notebook cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the plot backend to support interactive plotting\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import data and plot\n",
    "\n",
    "The following cell imports the NMR intensity matrix and chemical shift range $\\delta$ (*x-axis*), in part-per-million (ppm). We also import another matrix (Y), which contains metadata for each observation.\n",
    "\n",
    "### Metadata\n",
    "Y1 (first column of Y matrix) - represents the genotype (1: wild-type, 2: *sod-2* mutants, in original Y data matrix)\n",
    "\n",
    "Y2 (second column of Y matrix) - represents the age (1: younger L2 worms, 2: L4 worms, in original Y data matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "X = np.genfromtxt(\"./data/X_spectra.csv\", delimiter=',', dtype=None)\n",
    "Y = pds.read_csv(\"./data/worm_yvars.csv\",delimiter=',',dtype=None, header=None)\n",
    "ppm = np.genfromtxt(\"./data/ppm.csv\", delimiter=',', dtype=None)\n",
    "\n",
    "# Use pandas Categorical type to generate the dummy enconding of the Y vector (0 and 1) \n",
    "Y1 = pds.Categorical(Y.iloc[:, 0]).codes\n",
    "Y2 = pds.Categorical(Y.iloc[:, 1]).codes\n",
    "\n",
    "print('The X matrix contains {0} rows (observations/samples) and {1} columns (variables)'.format(X.shape[0], X.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: To apply the analyses exemplified in this notebook to any other dataset, just modify the cell above to import the data matrices and vectors X and Y from any other source file.\n",
    "\n",
    "The expected data types and formatting for **X** and **Y** are:\n",
    "\n",
    "   **X**: Any data matrix with n rows (observations/samples) and p columns (variables/features). The matrix should be provided as a [numpy.ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) object, with 2 dimensions, and with shape = (n, p). We recommend using the *numpy* function [numpy.genfromtxt](https://numpy.org/devdocs/reference/generated/numpy.genfromtxt.html) or the *pandas* [pandas.read_csv](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html) function to read the data from a text file. When using the *pandas.read_csv* function, extract the data matrix as a *numpy.ndarray* from the pandas.DataFrame object using the `.values` attribute. \n",
    "```\n",
    "X_DataFrame = pds.read_csv(\"./data/X_spectra.csv\")\n",
    "X = X_DataFrame.values\n",
    "```\n",
    "   \n",
    "   **Y** vectors: Each **Y** vector should be a 1-dimensional [numpy.ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) object, with a number and ordering of elements matching the rows in **X**. For continuous variables, any regular *numpy.ndarray* with a data type of `int` (integers only) or `float` can be used.\n",
    "   ```\n",
    "   Y_continuous = numpy.ndarray([23.4, 24, 0.3, -1.23], dtype='float')\n",
    "   ```\n",
    "To encode binary class labels, a *numpy.ndarray* of dtype `int`, with 0 and 1 as labels (e.g., 0 = Control, 1 = Case) must be used. The way in which classes are encoded will affect the model interpretation: the class labeled as 1 is used as the \"positive/case\" class by the *pyChemometrics* objects.\n",
    "   \n",
    "   In the example above, we used the *pandas* [Categorical](https://pandas.pydata.org/pandas-docs/stable/user_guide/categorical.html) datatype to handle the conversion of the original numerical values (1, 2) to the required (0, 1) labels. After converting a column to a `Categorical` datatype, the `.codes` attribute returns a vector with the same length of the original Y, but where each value is replaced by their integer (`int`) code. The correspondence between code and category can be inspected with the `categories` attribute. The order of the labels in `.codes` is the same as the order of the `categories` attribute (i.e. 0 is the first element in `categories`, 1 the second and so on).\n",
    "   ```\n",
    "   Y1 = pds.Categorical(Y.iloc[:, 1])\n",
    "   Y1.codes # The numerical label\n",
    "   Y1.categories # Original text or numerical description of the category\n",
    "   ```\n",
    "   [get_dummies](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.get_dummies.html) is another helpful function to perform dummy (0-1) encoding of variables. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the raw data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot in the next cell displays the raw NMR data. The X axis corresponds to the chemical shift in $\\delta_H$ ppm (part-per-million). The TSP chemical shift reference and the water peak region ($\\delta_H$ 5-4.6 ppm) have been previously removed, as these regions are not biologically informative. Since this is an NMR spectroscopy dataset we follow the convention of presenting the $\\delta$ppm scale with 0 on the right of the axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(ppm, X.T)\n",
    "plt.title(\"Raw data\")\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalisation\n",
    "\n",
    "The main purpose of normalisation is to account for systematic experimental differences in the overall intensity of observations (rows of the data matrix) that can confound the analysis. For example, different amount of material per sample, urinary dilution or an instrumental intensity drift caused by a batch or run-order effect. \n",
    "\n",
    "One of the first decisions in the statistical analysis of metabolic profile data is whether or not to apply \n",
    "normalisation, and which method to use. Mathematically, normalisation procedures estimate and apply an individual scaling factor to each observation in the data matrix (row-wise scaling of the data matrix).\n",
    "\n",
    "Here we will demonstrate and use the Probabilistic Quotient Normalisation ([Dieterle et al. 2006](https://www.ncbi.nlm.nih.gov/pubmed/16808434)) method. \n",
    "\n",
    "The PQN algorithm requires the selection of reference spectrum. A robust default choice is to select the median or mean spectrum as normalisation reference.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get reference median spectra\n",
    "median_spectra = np.median(X, axis=0)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Median Spectra\")\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.plot(ppm, median_spectra)\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After selecting a reference, each data point (variable) in each observation (row) is divided by the corresponding value of the same variable in the reference, to obtain a matrix of ratios to reference (*foldChangeMatrix*). The PQN coefficients (*pqn_coef*) corresponding to each observation are then estimated by taking the median across all variables row-wise. These coefficients are then used to normalised the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Calculate the ratio to the reference per variable for all spectrum, to obtain a matrix of \"fold changes\"\n",
    "foldChangeMatrix = X / median_spectra[None, :]\n",
    "\n",
    "# Estimate the median fold change per row (spectrum), which is the PQN normalisation coefficient. \n",
    "# We use the nanmedian function to automatically ignore NA's when estimating this factor\n",
    "pqn_coef = np.nanmedian(foldChangeMatrix, axis=1)\n",
    "\n",
    "# Normalise the spectra using the median ratio\n",
    "norm_X = X / pqn_coef[:, None]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the PQN coefficients, to assess the range of overall intensity/dilution/concentration ranges in the dataset which will be corrected by the normalisation procedure. \n",
    "\n",
    "In the next histogram, these have been re-expressed as a scaling factor, $\\frac{1}{PQN coefficient}$, for ease of interpretation: values smaller than one represent samples with an overall intensity higher than the median which will be downscaled, and the opposite for samples with value above 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(1/pqn_coef, bins=25)\n",
    "plt.xlabel(\"1/PQN Coefficient\")\n",
    "plt.ylabel('Counts')\n",
    "plt.title(\"Distribution of Normalisation factors\")\n",
    "plt.show()\n",
    "\n",
    "#Truncate extreme values to narrow histogram range\n",
    "sample_to_plot = 10\n",
    "idx_to_plot = ((foldChangeMatrix[sample_to_plot, :] <= 5) & (foldChangeMatrix[sample_to_plot, :] >= -5 ))\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Fold change to reference for sample: {0}'.format(sample_to_plot))\n",
    "plt.xlabel(\"Fold Change to median\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.hist(foldChangeMatrix[sample_to_plot, idx_to_plot], bins=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compare the original non-normalised data, plotted in blue, with the normalised spectra, in red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(ppm, X.T,  c='b', alpha=0.4)\n",
    "ax.plot(ppm, norm_X.T, c='r', alpha=0.4)\n",
    "ax.set_title(\"Normalised spectra\")\n",
    "ax.set_xlabel(\"$\\delta_H$ in ppm\")\n",
    "ax.set_ylabel(\"Intensity (a.u.)\")\n",
    "ax.invert_xaxis()\n",
    "ax.legend()\n",
    "plt.show()\n",
    "\n",
    "np.savetxt('./data/X_norm.csv', norm_X, delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scaling\n",
    "\n",
    "The choice of scaling is an important parameter choice when modelling the data with PCA and PLS regression algorithms. The scaling parameter impacts on the predictive performance of the model, the trends recovered from the data and model interpretation (i.e., loading vectors).\n",
    "\n",
    "Scaling consists in dividing each variable (column of the data matrix) by a constant value. Variables with a higher intensity tend to have a higher variance. Since PCA and PLS maximize the recovery of variance and covariance, on non-scaled datasets the model will be biased towards features with higher signal intensity. Biologically, it is more reasonable to ensure all variables should be be given equal weight independent of their signal intensity. Scaling the data matrix can remove this bias and enhance recovery of important trends from lower intensity signals.\n",
    "\n",
    "For more information about data scaling and transformation in metabolomics, we recommend [Berg et al. 2006](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1534033/)\n",
    "\n",
    "The following 3 types of scaling choices are commonly used with multivariate PCA and PLS analyses:\n",
    "\n",
    "Mean centring (MC) = $\\frac{X - \\mu}{\\sigma^{0}}$, *scale_power* = 0. Mean centre all variables and apply no scaling.\n",
    "\n",
    "Pareto scaling = $\\frac{X - \\mu}{\\sigma^{1/2}}$, *scale_power* = 1/2. Mean centre all variables and divide each variable by the square root of its standard deviation ($\\sigma$)\n",
    "\n",
    "Unit variance (UV) scaling = $\\frac{X - \\mu}{\\sigma^{1}}$, *scale_power* = 1. Mean center all variables and scale each variable by its standard deviation.\n",
    "\n",
    "\n",
    "In the *pyChemometrics.ChemometricsScaler* objects we use the *scale_power* argument to represent the exponent of standard deviation used in scaling.\n",
    "\n",
    "These three kinds of scaling are demonstrated below, followed by two other data transformation options, the logarithm and the square root transformation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) Mean-centring only (no scaling)\n",
    "\n",
    "Removal of the column-wise mean from the data matrix (each variable is recentered around its own mean). Mean centring is usually suggested (although not a hard-requirement, the data). \n",
    "\n",
    "We will apply the different scaling methods using the *pyChemometrics.ChemometricsScaler* object. We start by initialising a *ChemometricsScaler* object. The *fit* method extracts the scaling parameters (mean and standard deviation vectors) from a data matrix, while the *transform* method takes a data matrix as input and applies to it the scaling using the parameters estimated when the *fit* method was last used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc_scaler = ChemometricsScaler(scale_power=0)\n",
    "mc_scaler.fit(norm_X)\n",
    "mc_X = mc_scaler.transform(norm_X)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell we plot the mean centred data matrix. The main features of the spectrum are still relatively easy to pinpoint (which will aid interpretation of multivariate parameters), but the wide range of values still means that higher intensity variables will dominate the analysis. This *penalises* the detection of trends in low intensity signals, which is not, but also ensures that spurious variation in noise peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mc_X = preprocessing.scale(norm_X, axis=0, with_mean=True, with_std=False)\n",
    "plt.figure()\n",
    "plt.plot(ppm, mc_X.T)\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.title(\"Mean-centered scaling\")\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) Unit-variance (UV) scaling: Mean-centring followed by division of each variable (column) by its own standard deviation\n",
    "\n",
    "The spectroscopic peaks and signals of the spectrum are now barely recognisable (difficult interpretation), but now the latent variables detected will not be biased towards variation in the high intensity variables, which biologically is a more reasonable criteria. Low intensity noise values are more likely to be picked up, including baseline artefacts, and the magnitude of most variables is now equalised, making it harder to recognise the spectral profile. \n",
    "\n",
    "In the next example we use the *fit_transform* method, which learns the scaling parameters and returns a scaled data matrix in one go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uv_scaler = ChemometricsScaler(scale_power=1)\n",
    "uv_X = uv_scaler.fit_transform(norm_X)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Univariate and mean-centred\n",
    "#uv_X = preprocessing.scale(norm_X, axis=0, with_mean=True, with_std=True)\n",
    "plt.figure()\n",
    "plt.plot(ppm, uv_X.T)\n",
    "plt.title(\"Univariate scaling\")\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) Pareto scaling: Mean-centring followed by division of each variable (column) by the square root of its own standard deviation\n",
    "\n",
    "Pareto scaling provides a balance between the mean centering (no scaling) and UV scaling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "par_scaler = ChemometricsScaler(scale_power=1/2)\n",
    "pa_X = par_scaler.fit_transform(norm_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#std_dev_sq_root = np.sqrt(mc_X.std(axis=0))\n",
    "#pa_X = mc_X / std_dev_sq_root\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(ppm, pa_X.T)\n",
    "plt.title(\"Pareto scaling\")\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4) Logarithmic and square root transformation\n",
    "\n",
    "Another option is apply non-linear transformations to the data matrix, such as the logarithmic transform and the square-root transformation. \n",
    "\n",
    "Log-transformation is a non-linear transformation, which has the advantage of removing the effect of heteroskedasticity between and within variables, when the coefficient of variation of the methods is constant (variance increasing with the mean). The scaling effect it has on variables is somewhat similar to Pareto scaling. They are also recommended for univariate analysis. \n",
    "These transformations functions are not defined for 0 values. To perform logarithmic transform  we first add an offset to ensure all data points are non-negative and non-zero. Here it is calculated by adding the minimum value in the dataset plus 1. Any offset or other transformation should be carefully recorded. \n",
    "\n",
    "Logarithmic and square-root transformations are not implemented in the *pyChemometrics.ChemometricsScaler* object, but can be applied by transforming the data as shown below. After applying these transformations, the data should still be mean centred before PCA and PLS modeling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Log scaling\n",
    "\n",
    "# Offset is required to ensure all datapoints are non-negative and non-zero.\n",
    "# Here it is being calculated by simply adding 1 \n",
    "#to the minimum value in the dataset (to account for negative values in the noise region as well as 0).\n",
    "\n",
    "offset = np.min(norm_X) + 1\n",
    "log_X = np.log(norm_X + offset)\n",
    "mean_logvec = np.mean(log_X, axis=0)\n",
    "mclogX = (log_X - mean_logvec)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(ppm, mclogX.T)\n",
    "plt.title(\"Logarithmic transformation\")\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()\n",
    "\n",
    "# Square root transformation\n",
    "\n",
    "# Offset is required to ensure all datapoints are non-negative and non-zero.\n",
    "# Here it is being calculated by simply adding 1 \n",
    "#to the minimum value in the dataset (to account for negative values in the noise region as well as 0).\n",
    "\n",
    "offset = np.min(norm_X) + 1\n",
    "sqrt_X = np.sqrt(norm_X + offset)\n",
    "mean_sqrvec = np.mean(sqrt_X, axis=0)\n",
    "mcsqrt_X = (sqrt_X - mean_sqrvec)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(ppm, mcsqrt_X.T)\n",
    "plt.title(\"Square root transformation\")\n",
    "plt.xlabel(\"$\\delta_H$ in ppm\")\n",
    "plt.ylabel(\"Intensity (a.u.)\")\n",
    "plt.gca().invert_xaxis()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}