Preprocessing and Quality Control of LC-MS Data with the nPYc-Toolbox.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing and Quality Control of LC-MS Data with the nPYc-Toolbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial demonstrates how to use the LC-MS data processing modules of the nPYc-Toolbox, to import and perform some basic preprocessing and quality control of LC-MS data, and to output a final high quality dataset ready for data modelling. It is based on the quality control criteria previously described in [Lewis et al. 2016](https://www.ncbi.nlm.nih.gov/pubmed/27479709).\n",
    "\n",
    "Details of how to install all of the required dependencies and to set up your computing environment can be found here [Installing the nPYc-Toolbox](https://npyc-toolbox.readthedocs.io/en/latest/tutorial.html#Installing-the-nPYc-Toolbox).\n",
    "\n",
    "The dataset used in this example (DEVSET U RPOS xcms) is comprised of six samples of pooled human urine, aliquoted, and independently prepared and measured by ultra-performance liquid chromatography coupled to reversed-phase positive ionisation mode spectrometry (LC-MS, RPOS). Each source sample was separately prepared and assayed thirteen times. A pooled QC sample (study reference, SR) and independent external reference (long-term reference, LTR) of a comparable matrix was also acquired to assist in assessing analytical precision. See [Tutorial Datasets](https://npyc-toolbox.readthedocs.io/en/latest/tutorial.html#tutorial-datasets) and the Metabolights Study [MTBLS694](https://www.ebi.ac.uk/metabolights/MTBLS694) for details of the study, and [Recommended Study Design Elements](https://npyc-toolbox.readthedocs.io/en/latest/studydesign.html) for details of the various QC samples acquired.\n",
    "\n",
    "The [nPYc-toolbox-tutorials](https://github.com/phenomecentre/nPYc-toolbox-tutorials) contain all of the data required to run the tutorial Jupyter notebooks, full details of which are given [here](https://npyc-toolbox.readthedocs.io/en/latest/tutorial.html#preprocessing-and-quality-control-of-lc-ms-data-with-the-npyc-toolbox)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Import the nPYc-Toolbox and Configure the Jupyter Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Import the nPYc-Toolbox\n",
    "import nPYc\n",
    "\n",
    "# Import enumerations for sample type\n",
    "from nPYc.enumerations import VariableType, DatasetLevel, AssayRole, SampleType\n",
    "\n",
    "# Import normalisation objects for data normalisation\n",
    "from nPYc.utilities.normalisation import NullNormaliser, TotalAreaNormaliser, ProbabilisticQuotientNormaliser\n",
    "\n",
    "# Import matplotlib plotting, configure the Jupyter notebook to plot inline\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# Set up plotly to work in offline mode with the notebook\n",
    "from plotly.offline import init_notebook_mode, iplot\n",
    "init_notebook_mode(connected=True)\n",
    "\n",
    "# Set up to hide warnings (particularly Depreciation, RunTime warnings, these are not for the user to worry about!)\n",
    "# These lines can be commented to show warnings if of interest to advanced users\n",
    "import warnings\n",
    "warnings.filterwarnings(action='ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Import LC-MS Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to import the LC-MS data into an nPYc-Toolbox [Dataset](https://npyc-toolbox.readthedocs.io/en/latest/objects.html) object.\n",
    "\n",
    "The nPYc-Toolbox works on feature extracted LC-MS data, where the raw analytical data has been peak-picked into a one dimensional list of detected features, each of which is characterised by abundance and observed m/z and retention time. There are a number of peak-detection algorithms ([Spicer et al. 2017](https://www.ncbi.nlm.nih.gov/pubmed/28890673)), many of which are supported by the nPYc-Toolbox, however in this example, untargeted feature extraction was carried out using [XCMS](https://bioconductor.org/packages/release/bioc/html/xcms.html) [(Smith et al. 2006)](https://www.ncbi.nlm.nih.gov/pubmed/16448051).\n",
    "\n",
    "In brief, the raw data files were first converted from .RAW format to .mzML format using the msConvert tool from the ProteoWizard tookit ([Chambers et al. 2012](https://www.ncbi.nlm.nih.gov/pubmed/23051804)), for feature-detection with XCMS. msconvert was configured to discard scans below an intensity threshold of 100 (--filter 'threshold absolute 100 most-intense'), and only retain scans from the first acquisition function (--filter 'scanEvent 1'). Subsequent peak detection by XCMS used the centwave algorithm configured with a noise threshold of 600, mass accuracy of 25ppm, and minimum and maximum peakwidth of 1.5 and 5 seconds. Feature grouping was performed using the 'density' method.\n",
    "\n",
    "The LC-MS XCMS DevSet Dataset is located in 'DEVSET U RPOS xcms.csv', preceding the file path (in quotes) with the letter r ensures that the path is interpreted exactly as written.\n",
    "\n",
    "The following line creates an object representing the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msData = nPYc.MSDataset(r'DEVSET U RPOS xcms.csv', fileType='XCMS', sop='GenericMS')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The “sop” parameter points to a configuration file (encoded in JSON format) which contains a set of parameters to use during data import and pre-processing, see [Configuration Files](https://npyc-toolbox.readthedocs.io/en/latest/configuration/configuration.html) for full details.\n",
    "\n",
    "The default configuration file for LC-MS data, 'GenericMS', contains the recommended parameters for import and quality control of both human urine and plasma/serum data. For a list of all the parameters for MS data, see Table 5 in [Built-in Configuration SOPs](https://npyc-toolbox.readthedocs.io/en/latest/configuration/builtinSOPs.html).\n",
    "\n",
    "If required, users can create new configuration files, or indeed amend the existing documents with their own values, however, any of the parameters present in these files can also be overwritten by passing values into the data import command directly, without having to modify or generate the configuration files themselves.\n",
    "\n",
    "For LC-MS data, the majority of the parameters relate to preprocessing at the feature filtering stage, and therefore discussed in more detail below, but as an example, to change the threshold for filtering features based on residual standard deviation (RSD), the argument \"rsdThreshold\" can be overridden in the following manner:\n",
    "\n",
    "```\n",
    "msData = nPYc.MSDataset('DEVSET U RPOS xcms.csv', fileType='XCMS', sop='GenericMS', rsdThreshold=20)\n",
    "```\n",
    "\n",
    "Each nPYc-Toolbox Dataset object contains a name that can be changed as shown in the next cell. This name will be used in the summary and visualisation reports and in the file names of the exported data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msData.name = 'nPYc LC-MS Tutorial dataset'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Import Sample Metadata and Match to Acquired Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "The default way to add sample metadata is to prepare a CSV file which follows the set of conventions as described in [CSV Template for Metadata Import](https://npyc-toolbox.readthedocs.io/en/latest/samplemetadata.html#csv-template-for-metadata-import) and match it with the acquired data using the 'addSampleInfo' method.\n",
    "\n",
    "Based on the corresponding entries in the sample metadata CSV file, the acquired samples are categorised into different types, where 'Study Samples (SS)' comprise the main core of the study, and the others are acquired for specific roles in characterising data quality. The main QC samples here are the 'Study Reference (SR)' samples, which comprise a pool of study samples and are used to assess analytical stability across the run, and the 'Long Term Reference (LTR)' samples, a QC sample external to the study. Additionally, the 'Serial Dilution (SRD)' samples are key for assessing linearity of response i.e. that features are measured accurately with respect to their true abundance. For interest we have also included a blank, for full details see [Recommended Study Design Elements](https://npyc-toolbox.readthedocs.io/en/latest/studydesign.html).\n",
    "\n",
    "Although optional, this is recommended in order to make optimal use of the quality control features and visualisations provided by the nPYc-Toolbox.\n",
    "\n",
    "An example CSV file is provided, as given in 'DEVSET U RPOS Basic CSV.csv':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msData.addSampleInfo(descriptionFormat='Basic CSV', filePath=r'DEVSET U RPOS Basic CSV.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in [Datasets](https://npyc-toolbox.readthedocs.io/en/latest/objects.html), the spectral data, sample metadata and feature metadata can be inspected directly using:\n",
    "\n",
    "```\n",
    "dataset.intensityData\n",
    "dataset.sampleMetadata\n",
    "dataset.featureMetadata\n",
    "\n",
    "```\n",
    "\n",
    "In addition to the 'Basic CSV' format, it is also possible to extract analytical parameters directly from the raw data files. In this example, we have restricted data input to the sample metadata CSV file owing to the large file sizes of the raw LC-MS data files. However, if you would like to use this method for the tutorial dataset, the raw files can be downloaded directly from [MTBLS694](https://www.ebi.ac.uk/metabolights/MTBLS694) and the following used to match and extract the analytical parameters:\n",
    "\n",
    "```\n",
    "msData.addSampleInfo(descriptionFormat='Raw Data', filePath=r'/path to raw data files')\n",
    "```\n",
    "\n",
    "See [Analytical Parameter Extraction](https://npyc-toolbox.readthedocs.io/en/latest/samplemetadata.html#analytical-parameter-extraction) for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Generate Quality Control Reports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nPYc-Toolbox offers a series of reports, pre-set visualisations comprised of text, figures and tables to describe and summarise the characteristics of the dataset, and help the user assess the overall impact of quality control decisions (i.e. excluding samples or features and changing filtering criteria). \n",
    "\n",
    "For full details see [Reports](https://npyc-toolbox.readthedocs.io/en/latest/reports.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Sample Summary Report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first report can be used to check the expected samples against those acquired, in terms of numbers, sample type, and any samples either missing from acquisition or not recorded in the sample metadata CSV file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msData, 'sample summary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, for this dataset there are five samples with entries in the sample metadata CSV file, but missing from acquisition. This allows the user to quickly assess the completeness of the dataset and, for example, investigate why these samples were missing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Summary Report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The feature summary report provides visualisations summarising the quality of the dataset with regards to quality control criteria previously described in [Lewis et al. 2016](https://www.ncbi.nlm.nih.gov/pubmed/27479709). These include both assessment of potential run-order and batch effects, and metrics by which feature quality can be assessed.\n",
    "\n",
    "In order, these consist of:\n",
    "- Feature abundances (Figure 1)\n",
    "- Sum of total ion count, TIC (Figures 2 and 3)\n",
    "- Correlation to dilution (Figures 4, 5 and 7)\n",
    "- Residual standard deviation, RSD (Figures 5, 7 and 9)\n",
    "- Chromatographic peak width, if available (Figure 8)\n",
    "- Ion map (Figure 10)\n",
    "\n",
    "For several of these parameters (for example, correlation to dilution, RSD), acceptable default values are pre-defined in the configuration SOP, see [Built-in Configuration SOPs](https://npyc-toolbox.readthedocs.io/en/latest/configuration/builtinSOPs.html). If different values are required, these can be set by the user in the SOP directly, or by updating the 'Attribute', either at import, or by subsequent direct modification in the pipeline (see [Datasets](https://npyc-toolbox.readthedocs.io/en/latest/objects.html) for more details).\n",
    "\n",
    "For full details of each of the above, see [Feature Summary Report: LC-MS Datasets](https://npyc-toolbox.readthedocs.io/en/latest/reports.html#feature-summary-report-lc-ms-datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msData, 'feature summary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, for this dataset, the features are of reasonable intensity (Figure 1), and show the expected distributions in both correlation to dilution (Figure 4, positively skewed) and RSD (Figure 6, negatively skewed). Although there is a slight decrease in overall feature intensity over the course of the run (Figure 2) this is normal for LC-MS datasets, it is not associated with detector voltage (Figure 3), correction of this (by Batch and Run-Order Correction) will be addressed below. From the ion map (Figure 10) a large band of features can be seen at the very early retention times, these are in the solvent front and should be excluded (see below)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exporting Reports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default all reports are output directly to the notebook (as above), however, if html copies are required these can be automatically saved to the save directory by adding the optional input argument \"destinationPath\".\n",
    "\n",
    "For example, to save to the path defined in \"saveDir\":  \n",
    "\n",
    "```\n",
    "saveDir = '/path to save outputs'\n",
    "nPYc.reports.generateReport(msData, 'feature summary', destinationPath=saveDir)\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Assess Batch and Run-Order Effects and Apply Correction if Necessary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nPYc-Toolbox provides a module specifically for detecting and correcting for per-feature run-order and batch effects in LC-MS datasets, by characterising the effect in reference samples and interpolating a correction factor to the intermediate samples. For full details see [Batch & Run-Order Correction](https://npyc-toolbox.readthedocs.io/en/latest/batchAndROCorrection.html).\n",
    "\n",
    "In brief, for each feature, a LOWESS estimator is fitted on a series of consecutive study reference samples (with the default number of samples, \"window=11\") for each analytical batch (which can be defined by the user), and this fit used to correct all samples in the run (excluding SRD samples). Batch differences are corrected by aligning median feature intensities in the study reference samples between batches.\n",
    "\n",
    "Initially, batch and run-order correction performance can be assessed on a subset of features prior to running on the whole dataset using the 'batch correction assessment' report:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msData, 'batch correction assessment')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By comparing the results across these surveyed features, the parameters for and necessity of correction can be assessed. As this is a critical step in preprocessing, here we give a full worked example of what to look for, there are four main questions which need to be addressed:\n",
    "\n",
    "**1. Is the window of the LOWESS smoother appropriate?**\n",
    "\n",
    "It is crucial that only broad and not narrow trends are being fitted, i.e. that the data is not being over-fitted. This can firstly be visualised by observing the fit of the LOWESS (yellow line) in the plots above, and secondly, by looking at the results post-correction using the 'batch correction summary' report.\n",
    "\n",
    "The default window over which the LOWESS is fitted is 11 study reference (SR) samples. In the following, to illustrate the effects of over-fitting, we can generate a corrected dataset with \"window=3\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "msDatacorrected = nPYc.batchAndROCorrection.correctMSdataset(msData, window=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To assess the fit and performance of correction, the 'batch correction summary' report can then be generated using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msData, 'batch correction summary', msDataCorrected=msDatacorrected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In terms of assessing over-fitting, the key figures here are Figures 2 and 4, in both of these, by comparing the pre-correction and post-correction plots it is clear that the relationship between the study reference and long-term reference samples has been completely skewed (in this case the RSD in SR samples post-correction (Figure 3) also cannot be calculated owing to the presence of zero values). The aim of correcting the data is to choose a window where any broad trends are corrected, but the local relationship between different sample types is maintained."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2. Does the correction need to be applied in different sub-batches?**\n",
    "\n",
    "If the samples have been acquired in one continuous run, it is likely that correction in one overall batch is appropriate. However, if there have been any substantial stoppages, or the samples have been run in multiple batches it may be appropriate to split the correction into a number of batches. Practically, this means that a LOWESS smoother is fitted to the study reference samples in each batch, then the batches are aligned to each other at the end using the median intensity of the reference samples. This prevents errors in correction which might otherwise occur around the start of a new batch or after a stoppage, owing to a step, rather than a gradual change in overall intensity.\n",
    "\n",
    "As this dataset contains only a small number of samples, and with no long stoppages, running in one batch would be suitable. However, if splitting the dataset into multiple batches were required, there are a number of helpful functions, described here:\n",
    "\n",
    "An interactive plot can first be generated to assess at which point to start a new batch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = nPYc.plotting.plotTICinteractive(msData, plottype='Sample Type')\n",
    "iplot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New batches can be added using 'amendBatches', where the input argument is the run order of the first sample required to be in the new batch (as determined from the interactive plot above).\n",
    "\n",
    "For example, in this case a new batch could be started at '136' using:\n",
    "\n",
    "```\n",
    "msData.amendBatches(136)\n",
    "\n",
    "```\n",
    "\n",
    "Re-generating the interactive plot would then show the study samples coloured in two different batches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3. Do any SR samples need to be excluded?** \n",
    "\n",
    "Generally speaking, if the window of the LOWESS smoother is set to an appropriate value, a single SR sample with an outlying TIC should not affect the results of the fitting.\n",
    "\n",
    "In this example dataset, there are five SR samples at the beginning, and five after the short break (from 'Run Order'=136) which have the same composition as all the other SR samples in the run, but were injected from a vial rather than from a well plate. In this case, if a different response was seen in the 'batch correction assessment' report, it could be appropriate to mask and permanently exclude these specific samples using:\n",
    "\n",
    "```\n",
    "msData.excludeSamples([58, 59, 60, 61, 62, 136, 137, 138, 139, 140], on='Run Order', message='Vial SR')\n",
    "msData.applyMasks()\n",
    "```\n",
    "\n",
    "For more information on these functions see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4. Is correction required?** \n",
    "\n",
    "If there is no observable trend in the run-order, and samples were acquired in one continuous batch it could be that correction is not required!\n",
    "\n",
    "The analytical 'Multivariate Report' can also be useful when making this decision (see below for details and an example).\n",
    "\n",
    "In this case, correction is clearly required (see the 'Multivariate Report' run on data before correction), and can be run with the default parameters using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "msDatacorrected = nPYc.batchAndROCorrection.correctMSdataset(msData)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, the 'batch correction summary' report can then be generated using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msData, 'batch correction summary', msDataCorrected=msDatacorrected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the default parameters it can be seen that the relationship between the SR and LTR samples is conserved, but the gradual decline in overall intensity during the run is corrected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Exclude Samples and/or Features if Required"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset objects contain two internal 'mask' vectors, the 'sampleMask' and 'featureMask', which store whether a sample or feature respectively should be used when calculating QC metrics, visualised in the reports and finally exported, see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html).\n",
    "\n",
    "There are several functions which modify these masks, which are useful at various stages of quality control and in preparing a final dataset for export."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mask Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For LC-MS datasets, features should be filtered based on their individual precision and accuracy to ensure that only high quality features are retained. \n",
    "\n",
    "The nPYc-Toolbox has been developed based on the quality control criteria previously described in [Lewis et al. 2016](https://www.ncbi.nlm.nih.gov/pubmed/27479709), and uses three key criteria for feature filtering, see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html):\n",
    "\n",
    "1. Correlation to dilution (calculated from SRD samples), assess that features are measured accurately with respect to their true abundance\n",
    "2. RSD (calculated from SR samples), assesses that features are measured precisely from multiple acquisitions of the same sample across the run\n",
    "3. Comparing RSD in SR and SS, the variance in SS should always be greater than that seen in the SR samples\n",
    "\n",
    "The distribution of correlation to dilution, and RSD can be visualised using the 'Feature Summary Report'.\n",
    "\n",
    "A report summarising the number of features passing selection with different criteria can also be produced using: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'feature selection')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this report, the total number of features passing selection (with the default criteria) is 9920.\n",
    "\n",
    "If more (or less) stringent criteria are required these can be set in the object 'Attributes', for example:\n",
    "\n",
    "To amend the RSD threshold: \n",
    "\n",
    "```\n",
    "msDatacorrected.Attributes['rsdThreshold'] = 20\n",
    "```\n",
    "\n",
    "Or to amend the correlation to dilution threshold:\n",
    "\n",
    "```\n",
    "msDatacorrected.Attributes['corrThreshold'] = .8\n",
    "```\n",
    "\n",
    "For more details and a full list of parameters see [Built-in Configuration SOPs](https://npyc-toolbox.readthedocs.io/en/latest/configuration/builtinSOPs.html).\n",
    "\n",
    "Using the 'updateMasks' function, these criteria will be checked and the featureMask updated to mask those features not meeting the thresholds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "msDatacorrected.updateMasks(filterFeatures=True, filterSamples=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results of masking can be visualised using the 'feature summary' report.\n",
    "\n",
    "Using \"withExclusions=True\" means the report is generated as if any masked features were excluded from the dataset, which allows assessment of the results of filtering before the features are permanently excluded from the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'feature summary', withExclusions=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mask Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By setting preferences with the 'sampleTypes' and 'assayRoles' argument, we can limit our dataset to contain only samples of interest, masking all samples which do not fall into these categories.\n",
    "\n",
    "For example, to limit our dataset to contain only study samples ('SampleType.StudySample, AssayRole.Assay') and study reference samples ('SampleType.StudyPool, AssayRole.PrecisionReference') we would run the following (see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html) and [Enumerations](https://npyc-toolbox.readthedocs.io/en/latest/enumerations.html) for more details).\n",
    "\n",
    "As we have already set up the feature mask, we can set \"filterFeatures=False\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msDatacorrected.updateMasks(sampleTypes=[SampleType.StudySample, SampleType.StudyPool], assayRoles=[AssayRole.Assay, AssayRole.PrecisionReference], filterFeatures=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Long term reference, sample blank and dilution series samples have now been masked (marked for exclusion).\n",
    "\n",
    "**Important note**, here we have masked features prior to masking samples, the order is important in some cases, for example when filtering features by correlation to dilution the dataset must contain the serial dilution samples.\n",
    "\n",
    "The results of masking can be summarised using the 'sample summary' report."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'sample summary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excluding Specific Samples and/or Features "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 'updateMasks' function works to mask samples or features not meeting specific criteria, in addition to this, the nPYc-Toolbox also contains two additional methods to mask specific samples or features directly, 'excludeSamples' and 'excludeFeatures' respectively, see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html).\n",
    "\n",
    "Each of these funtions takes three input arguments; firstly, a list of sample or feature identifiers; secondly, the name of the column in 'sampleMetadata' (for 'excludeSamples') or 'featureMetadata' (for 'excludeFeatures') where these identifiers can be found; and finally an optional message as to why these samples or features have been flagged for exclusion.\n",
    "\n",
    "For example, to exclude a theoretical sample with an outlying TIC at 'Run Order' '93':\n",
    "\n",
    "```\n",
    "msDatacorrected.excludeSamples([93], on='Run Order', message='outlying TIC')\n",
    "```\n",
    "\n",
    "In the feature dimension, for LC-MS assays, there are regions of retention time outside of the useful range. For the reverse-phase positive assay this corresponds to features with a retention time of less than 0.6 minutes, or greater than 10.5 minutes (as seen from the 'feature summary' report ion map above). These can be masked using the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msDatacorrected.excludeFeatures(msDatacorrected.featureMetadata[msDatacorrected.featureMetadata['Retention Time'] < 0.6]['Feature Name'], on='Feature Name', message='Outside RT limits')\n",
    "msDatacorrected.excludeFeatures(msDatacorrected.featureMetadata[msDatacorrected.featureMetadata['Retention Time'] > 10.5]['Feature Name'], on='Feature Name', message='Outside RT limits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, the results of masking can be visualised using the 'feature summary' report, again with \"withExclusions=True\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'feature summary', withExclusions=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Permanently Exclude Masked Samples/Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once satisfied with the sample and feature masks, exclusions can be applied (permanently removed from the dataset) using the 'applyMasks' function.\n",
    "\n",
    "This method should be used only when it is absolutely certain that the masked features and samples are to be removed, as the excluded data will otherwise have to be re-imported.\n",
    "\n",
    "Before masks have been applied, however, feature/sample masking can be changed by first re-setting the masks to include all samples/features:\n",
    "\n",
    "```\n",
    "msDatacorrected.initialiseMasks() \n",
    "```\n",
    "\n",
    "Then different feature/sample exclusions can be applied as required.\n",
    "\n",
    "For details see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html)\n",
    "\n",
    "In this case we are happy with the masking, and so features and samples can be permanently excluded using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msDatacorrected.applyMasks()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Analytical Multivariate Quality Control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nPYc-Toolbox provides the capacity to generate a principal component analysis (PCA) model of the data (via the [pyChemometrics](https://github.com/phenomecentre/pyChemometrics) module), and subsequently, to use this to assess data quality, identify potential sample and feature outliers, and determine any potential analytical associations with the main sources of variance in the data ([Multivariate Analysis](https://npyc-toolbox.readthedocs.io/en/latest/multivariate.html))\n",
    "\n",
    "A PCA model can be generated using 'exploratoryAnalysisPCA', and there are a number of parameters which can be optimised depending on the dataset (see [PCA Model](https://npyc-toolbox.readthedocs.io/en/latest/multivariate.html#pca-model) for full details).\n",
    "\n",
    "One key parameter is 'scaling', which divides each column in the data matrix by its respective standard deviation raised to a power of the scaling parameter. This parameter can range in value between 0 and 1, and recommended values are 0 for mean-centering only, 0.5 for Pareto scaling and 1 for unit variance (UV) scaling. The outcome of PCA model will vary based on the scaling method selected, and different scaling functions can be appropriate depending on the data itself and the question being asked of the data, see [van der Berg et al. 2006](https://www.ncbi.nlm.nih.gov/pubmed/16762068)\n",
    "\n",
    "The default scaling is unit variance (\"scaling=1\"), which scales every variable to have a variance of one, and thus all variables (despite their different magnitudes) become equally important in the model.\n",
    "\n",
    "Each model is cross-validated using 7-fold cross-validation and the recommended number of principal components automatically estimated based on two criteria, when either one of these is met no more components will be added and the PCA model will be returned. There criteria are:\n",
    "1. \"minQ2\": Q2 is the variance predicted by each component (from cross-validation), when adding a component does not improve Q2 by at least this value (default \"minQ2=0.05\") then no more components will be added.\n",
    "2. \"maxComponents\": this defines the maximum number of components (default \"maxComponents=10\") returned by the model (regardless of Q2 increases).\n",
    "\n",
    "Again, these parameters can be amended by adding them as input arguments to 'exploratoryAnalysisPCA'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PCAmodel = nPYc.multivariate.exploratoryAnalysisPCA(msDatacorrected, scaling=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytical multivariate report provides visualisations summarising the largest sources of variance in the dataset (from the PCA model generated) with particular emphasis on any potential analytical sources, as defined in 'analyticalMeasurements' in the [Built-in Configuration SOPs](https://npyc-toolbox.readthedocs.io/en/latest/configuration/builtinSOPs.html).\n",
    "\n",
    "These consist of:\n",
    "- Scree plot of variance (Figure 1)\n",
    "- Scores plots coloured by sample type (Figure 2)\n",
    "- Strong sample outliers (Figure 3)\n",
    "- DmodX sample outliers (Figure 4)\n",
    "- Loadings plots (Figure 5)\n",
    "- Distribution of analytical parameters (Figure 6)\n",
    "- Heatmap of potential associations between analytical parameters and the main sources of variance (Figures 7 and 8)\n",
    "- Scores plots coloured by analytical parameters with potential association (Figures 9-11)\n",
    "\n",
    "For full details of each of the above, see [Multivariate Analysis Report](https://npyc-toolbox.readthedocs.io/en/latest/multivariate.html#multivariate-analysis-report)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nPYc.reports.multivariateReport(msDatacorrected, PCAmodel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the SR samples cluster tightly in the PCA scores plots, and no strong associations are observed in the heatmaps between analytical parameters and samples scores (main sources of variance in the data), we can conclude that the data is of high quality, and ready to be exported.\n",
    "\n",
    "For interest, we can also run the multivariate report on the dataset before correction, first filtering the samples to match that of the corrected dataset, and the features to keep only those passing the quality control criteria:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Filter samples (keeping SR and SS only):\n",
    "msData.updateMasks(sampleTypes=[SampleType.StudySample, SampleType.StudyPool], assayRoles=[AssayRole.Assay, AssayRole.PrecisionReference], filterFeatures=True)\n",
    "\n",
    "# Generate PCA model (with exclusions set to True):\n",
    "PCAmodel_precorrection = nPYc.multivariate.exploratoryAnalysisPCA(msData, scaling=1, withExclusions=True)\n",
    "\n",
    "# Generate report:\n",
    "nPYc.reports.multivariateReport(msData, PCAmodel_precorrection, withExclusions=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, from the heatmap of correlation to analytical parameters (Figure 7), and the corresponding analytical parameter coloured scores plots (Figure 9) it can be seen that run order correction is really necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scores and loadings plots can also be explored interactively with the 'plotScoresInteractive' and 'plotLoadingsInteractive' functions.\n",
    "\n",
    "**Interactive scores plot**\n",
    "\n",
    "For example, to plot the scores plot for principal component 1 vs. principal component  2 (\"components=[1, 2]\") with points coloured by values in msDatacorrected.sampleMetadata['Class'] (the colour is defined by the third input argument and can be any column name in the sample metadata):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = nPYc.plotting.plotScoresInteractive(msDatacorrected, PCAmodel, 'Class', components=[1, 2])\n",
    "iplot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, 'Class' gives information on which sample the spectra relate to (see 'DEVSET U RPOS Basic CSV.csv' and [Tutorial Datasets](https://npyc-toolbox.readthedocs.io/en/latest/tutorial.html#tutorial-datasets) for full details). Spectra from multiple acquisitions of the same sample cluster strongly, with the 'Study Pool' (Study Reference QC samples) clustering in the centre of the plot and clear separation can be seen between the different samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Interactive loadings plot**\n",
    "\n",
    "Similarly, to plot the loadings, here for principal component 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = nPYc.plotting.plotLoadingsInteractive(msDatacorrected, PCAmodel, component=2)\n",
    "iplot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Finalise and Export Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once no further exclusions or preprocessing is required, the final dataset can be exported.\n",
    "\n",
    "The 'final report' compiles information about the samples acquired, and the overall quality of the dataset taking key figures from feature and multivariate reports to give an overall summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'final report', pcaModel=PCAmodel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subsequently, datasets can be exported in a variety of formats with the 'exportDataset' method (see [Exporting Data](https://npyc-toolbox.readthedocs.io/en/latest/exportingdata.html)).\n",
    "\n",
    "By default datasets are exported to the current working directory, however, if files are required to be exported to a defined path, this can be done by adding the optional input argument \"destinationPath\".\n",
    "\n",
    "For example, to save to the path defined in \"saveDir\":  \n",
    "\n",
    "```\n",
    "saveDir = '/path to save outputs'\n",
    "msDatacorrected.exportDataset(saveFormat='UnifiedCSV', destinationPath=saveDir)\n",
    "\n",
    "```\n",
    "\n",
    "To export a single CSV file, which contains a row for every sample, and a column for every feature, alongside all of the sample and feature specific metadata, set \"saveFormat=UnifiedCSV\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msDatacorrected.exportDataset(saveFormat='UnifiedCSV', destinationPath='.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, to export the 'final report':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'final report', pcaModel=PCAmodel, destinationPath='.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. OPTIONAL: Normalise Data and Export"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dilution effects on global sample intensity can be normalised by using one of the available [Normalisation classes](https://npyc-toolbox.readthedocs.io/en/latest/normalisation.html).\n",
    "\n",
    "This is especially critical in urine samples when samples can vary substantially in concentration. Here we apply Probabilistic Quotient Normalisation (PQN) [Deterle et al. 2006](https://pubs.acs.org/doi/10.1021/ac051632c).\n",
    "\n",
    "Setting the 'Normalisation' parameter in msDataCorrected automatically applies the selected normalisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msDatacorrected.Normalisation = ProbabilisticQuotientNormaliser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the normalised data exported (here we first change the name for clarity and so as not to overwrite existing data):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Append name\n",
    "msDatacorrected.name = msDatacorrected.name + '_PQN_normalised'\n",
    "\n",
    "# Export dataset\n",
    "msDatacorrected.exportDataset(destinationPath='.')\n",
    "\n",
    "# Export final summary report\n",
    "nPYc.reports.generateReport(msDatacorrected, 'final report', destinationPath='.')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}