Preprocessing and Quality Control of Targeted LC-MS Data (PeakPantheR) 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": [
    "*peakPantheR* is a software package implemented in R for extracting signals or regions of interest of known identification from raw LC-MS data, to provide a fully annotated dataset. It is available to download from [Bioconductor](https://bioconductor.org/packages/peakPantheR/) or [GitHub](https://github.com/phenomecentre/peakPantheR), where tutorial vignettes can also be found. \n",
    "\n",
    "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 *peakPantheR* generated 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) 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, details of which are given [here](https://npyc-toolbox.readthedocs.io/en/latest/tutorial.html)."
   ]
  },
  {
   "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 PeakPantheR Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to import the PeakPantheR data into an nPYc-Toolbox [Dataset](https://npyc-toolbox.readthedocs.io/en/latest/objects.html) object.\n",
    "\n",
    "In brief, the raw LC-MS 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-annotation with *peakPantheR*. 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'). *peakPantheR* was then used to extract the intensities of 103 annotations with associated regions of interest (ROI), windows of 'Retention Time' and 'm/z'. Once the ROIs for each annotation had been user optimised, the final data tables were generated by running the 'outputAnnotationResult' method.\n",
    "\n",
    "Running *peakPantheRs* 'outputAnnotationResult' produces a number of files including the following, which are required by the nPYc-Toolbox:\n",
    "- 'peakArea.csv' : integrated areas for each sample (row) and feature (column)\n",
    "- 'cpdMetadata.csv' : associated feature metadata, including compound name, retention time and m/z\n",
    "\n",
    "The following line creates an object representing the *peakPantheR* dataset, preceding the file path (in quotes) with the letter r ensures that the path is interpreted exactly as written."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msData = nPYc.MSDataset(r'DevSet_peakArea.csv', 'csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 PeakPantheR 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": [
    "For feature metadata, information is added from the 'cpdMetadata.csv' file using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msData.addFeatureInfo('DevSet_cpdMetadata.csv', featureId='cpdID')"
   ]
  },
  {
   "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",
    "msData.intensityData\n",
    "msData.sampleMetadata\n",
    "msData.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 16 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. In this case, for two study samples ('DevSet5_Rep13' and 'DevSet6_Rep13') and two QC samples ('External Reference 8' and 'Study Pool 8', initial data conversion from .RAW to .mzML failed, and the 12 sample 'Blanks' were user excluded for simplicity as they contained no signals for the 103 *peakPantheR* annotations."
   ]
  },
  {
   "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) the Retention Time and m/z distribution of the dataset can be observed."
   ]
  },
  {
   "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": false
   },
   "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, for a full worked example please see 'Preprocessing and Quality Control of LC-MS Data with the nPYc-Toolbox.ipynb'.\n",
    "\n",
    "Here, correction is required, and can be run with the default parameters using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msDatacorrected = nPYc.batchAndROCorrection.correctMSdataset(msData)"
   ]
  },
  {
   "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": [
    "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 99.\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",
    "```"
   ]
  },
  {
   "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 adding violin plots for the distribution by sample type for each annotated feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(msDatacorrected, 'final report peakpanther')"
   ]
  },
  {
   "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(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 peakpanther', 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 peakpanther', destinationPath='.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}