Import, scaling and transfomation.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 different types of scaling to the data matrix and visualise their impact on the raw data matrix and interpretation.\n",
    "\n",
    "<br>\n",
    "\n",
    "*Credits: This tutorial was originally created by Gonçalo Correia and was adapted by Lukas Kopecky and Frederico Soares in November 2023.*"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Preliminary Steps"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "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",
    "import pandas as pd\n",
    "from sklearn import preprocessing\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "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,
   "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"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "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",
    "We will now import the LC-MS data with the metadata (Y variables) and feature annotation for LC-MS.\n",
    "\n",
    "Then we split the data into two parts:\n",
    "\n",
    "**rpos_x_matrix** - LC-MS data matrix\n",
    "\n",
    "**gender_y** - Metadata that will act as the response variable for the PLS-DA model (sex/gender in this instance)\n",
    "\n",
    "Then we extract the feature annotations:\n",
    "\n",
    "**retention_times**, **mz_values** - annotation for the features of the rpos_x_matrix data\n",
    "\n",
    "<br>\n",
    "\n",
    "*NB - The dataset used in this tutorial has reduced features to speed up the model fitting proces. Full data available from [https://zenodo.org/doi/10.5281/zenodo.4053166](https://zenodo.org/doi/10.5281/zenodo.4053166). We recommend you to try to run the tutorial with the full dataset in your own time, or you can even try to run it using your own data.*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Load the dataset\n",
    "dementia_rpos_dataset = pd.read_csv(\"./data/Dementia U RPOS_combinedData.csv\",delimiter=',')\n",
    "\n",
    "# Inspect the dataset\n",
    "dementia_rpos_dataset.head(10)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Delete samples where outcome variable is unknown (In this example QC samples)\n",
    "dementia_rpos_dataset = dementia_rpos_dataset[~dementia_rpos_dataset['Gender'].isnull()]\n",
    "\n",
    "rpos_x_matrix = dementia_rpos_dataset.iloc[:, 5::].values\n",
    "\n",
    "X = dementia_rpos_dataset.iloc[:, 5:505].values  #Subset for demonstration \n",
    "\n",
    "variable_names = dementia_rpos_dataset.columns[5::]\n",
    "\n",
    "# Use pandas Categorical type\n",
    "Y = pd.Categorical(dementia_rpos_dataset['Gender']).codes\n",
    "Y2 = dementia_rpos_dataset['Age'].values\n",
    "\n",
    "# Extract the retention times and m/z to use in 2D plots of the dataset\n",
    "retention_times = np.array([x.split('_')[0] for x in variable_names], dtype='float')/60\n",
    "mz_values = np.array([x.split('_')[1][0:-3] for x in variable_names], dtype='float')"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "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 a series of histograms of distributions of selected variables. The X axis corresponds to the peak relative abundances amd the Y axis their frequencies.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(X).plot.hist(bins=100, legend=False)"
   ]
  },
  {
   "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(X)\n",
    "mc_X = mc_scaler.transform(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": [
    "pd.DataFrame(mc_X).plot.hist(bins=100, legend=False)"
   ]
  },
  {
   "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",
    "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(X)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Univariate and mean-centred\n",
    "pd.DataFrame(uv_X).plot.hist(bins=100, legend=False)"
   ]
  },
  {
   "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(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",
    "pd.DataFrame(pa_X).plot.hist(bins=100, legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 4) Logarithmic and square root transformation\n",
    "\n",
    "Another option is to apply non-linear transformations to the data matrix, such as the logarithmic transform and the square-root transformation. "
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Log1p scaling\n",
    "log1p_x = np.log1p(X)\n",
    "\n",
    "pd.DataFrame(log1p_x).plot.hist(bins=100, legend=False)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "\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. "
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "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(X) + 1\n",
    "log_X = np.log(X + offset)\n",
    "mean_logvec = np.mean(log_X, axis=0)\n",
    "mclogX = (log_X - mean_logvec)\n",
    "\n",
    "pd.DataFrame(mclogX).plot.hist(bins=100, legend=False)\n",
    "\n",
    "\n",
    "# The code above is equivalent to: \n",
    "# mc_scaler = ChemometricsScaler(scale_power=0)\n",
    "# mclogX = mc_scaler.fit_transform(np.log1p(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "\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(X) + 1\n",
    "sqrt_X = np.sqrt(X + offset)\n",
    "mean_sqrvec = np.mean(sqrt_X, axis=0)\n",
    "mcsqrt_X = (sqrt_X - mean_sqrvec)\n",
    "\n",
    "pd.DataFrame(mcsqrt_X).plot.hist(bins=100, legend=False)\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "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
}