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

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing and Quality Control of NMR Data with the nPYc-Toolbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial demonstrates how to use the NMR data processing modules of the nPYc-Toolbox, to import and perform some basic preprocessing and quality control of NMR data, and to output a final high quality dataset ready for data modelling. It is based on the quality control criteria previously described in [Dona et al. 2014](https://www.ncbi.nlm.nih.gov/pubmed/25180432).\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 1D NMR) is comprised of six samples of pooled human urine, aliquoted, and independently prepared and measured by 1H NMR spectroscopy. 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-nmr-targeted-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 and Preprocess NMR Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second step is to import the 1D NMR raw data files (Bruker format) into a nPYc-Toolbox [Dataset](https://npyc-toolbox.readthedocs.io/en/latest/objects.html) object.\n",
    "\n",
    "The \"rawDataPath\" parameter sets the location of the NMR raw data.\n",
    "\n",
    "The NMR DevSet Dataset is located in 'DEVSET U 1D NMR raw data files'. This folder contains 94 directories, each corresponding to a spectrum acquired with the ‘noesygppr1d’ pulse sequence. The Fourier Transform, apodization and phasing have already been performed with the vendor software (i.e. TopSpin). \n",
    "\n",
    "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": [
    "rawDataPath = r'./DEVSET U 1D NMR raw data files'"
   ]
  },
  {
   "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 nPYc-Toolbox contains two default configuration files, 'GenericNMRurine' and 'GenericNMRblood'. These contain the recommended parameters for import and quality control of human urine and plasma/serum biofluid, respectively. For a list of all the parameters for NMR data, see the 'NMRDataset Objects' table in [Built-in Configuration SOPs](https://npyc-toolbox.readthedocs.io/en/latest/configuration/builtinSOPs.html).\n",
    "\n",
    "Since this is a urine biofluid dataset, we will use the 'GenericNMRUrine' file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sop = 'GenericNMRurine'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"pulseProgram\" parameter defines the specific NMR experiment pulse program to import, in this case 'noesygppr1d' - a standard 1D experiment with a NOE water pre-saturation. \n",
    "\n",
    "The specific text set in \"pulseProgram\" depends on the name of the pulse program (PULPROG) set when acquiring the data, and should match this exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pulseProgram = 'noesygppr1d'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following line creates an object representing the dataset and triggers pre-processing of the NMR spectra, including calibration to a reference peak, and interpolation of all spectra onto a common scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nmrData = nPYc.NMRDataset(rawDataPath, pulseProgram=pulseProgram, sop=sop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If required, users can create new configuration files, or indeed amend the existing documents with their own values, see [Configuration Files](https://npyc-toolbox.readthedocs.io/en/latest/configuration/configuration.html), 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 example, to interpolate the spectra to a higher resolution than the default, the argument \"variableSize\" can be overridden in the following manner:\n",
    "\n",
    "```\n",
    "nmrData = nPYc.NMRDataset(rawDataPath, pulseProgram=pulseProgram, sop=sop, variableSize=64000)\n",
    "```\n",
    "\n",
    "Each nPYc 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": [
    "nmrData.name = 'nPYc NMR Tutorial dataset'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Import Sample Metadata and Match to Acquired Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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' 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' samples, which comprise a pool of study samples and are used to assess analytical stability across the run. For interest we have also included some 'Long Term Reference' samples (a QC sample external to the study) and 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 1D NMR Basic CSV.csv:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nmrData.addSampleInfo(descriptionFormat='Basic CSV', filePath=r'DEVSET U 1D NMR 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",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "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": {},
   "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(nmrData, 'sample summary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, for this dataset there are two 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.\n",
    "\n",
    "From the 'sample summary' report, it can be seen that corresponding information is missing from the sample metadata CSV file for one sample. This sample is listed as having unknown type, and missing information. As above, this allows the user to quickly determine whether information should be added to the sample metadata CSV file for this sample, or whether the spectrum should be excluded from the final dataset (see below)."
   ]
  },
  {
   "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 [Dona et al. 2014](https://www.ncbi.nlm.nih.gov/pubmed/25180432).\n",
    "\n",
    "In order, these consist of:\n",
    "- Chemical shift calibration (Figure 1)\n",
    "- Line width (Figures 2)\n",
    "- Baseline consistency (Figure 3)\n",
    "- Quality of solvent suppression (Figure 3)\n",
    "\n",
    "For each parameter, 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 (as above), or by subsequent direct modification in the pipeline (see [Datasets](https://npyc-toolbox.readthedocs.io/en/latest/objects.html) for more details).\n",
    "\n",
    "Any samples failing any of the above criteria are flagged at the end of the report.\n",
    "\n",
    "For full details of each of the above, see [Feature Summary Report: NMR Datasets](https://npyc-toolbox.readthedocs.io/en/latest/reports.html#feature-summary-report-nmr-datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(nmrData, 'feature summary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, for this dataset, there is only one sample which fails the quality control criteria on water suppression quality.\n",
    "\n",
    "Using the interactive plotting function \"plotSpectraInteractive\" we can plot this spectrum to help decide whether it should indeed be excluded (and potentially re-acquired), or whether the solvent suppresion region for exclusion could be extended and the sample kept in the dataset (see below for an example).\n",
    "\n",
    "Note. \"plotSpectraInteractive\" works best with a small number of spectra. The specific samples to plot can be specified by the \"samples\" argument, here set to plot the first 5 spectra (with indices 0-4) and the sample failing on solvent peak suppression (with index 53 - index given in Table 1 at the end of the feature summary report above). The \"sampleLabels\" argument allows the user to set the labelling for the spectra, here we use 'Sample File Name', but 'Sample ID' could also be appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iplot(nPYc.plotting.plotSpectraInteractive(nmrData, samples=[0, 1, 2, 3, 4, 53], sampleLabels='Sample File Name'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By zooming in on the water suppression region, it can be seen that this spectrum has a significantly larger affected region. Depending on the potential value of peaks in this region, the solvent suppression region for exclusion could be extended, however, in order to get back to a normal baseline for this spectrum we would need to cut from 4-5 ppm, as such, in this case, we simply exclude this sample (see the next section)."
   ]
  },
  {
   "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(nmrData, 'feature summary', destinationPath=saveDir)\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. 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 Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 'updateMasks' function can be used to automatically mask samples (and/or features).\n",
    "\n",
    "For now, we do not mask any features, setting \"filterFeatures=False\".\n",
    "\n",
    "By default, any samples of unknown 'SampleType' and/or 'AssayRole' will be masked.\n",
    "\n",
    "In addition, samples failing specific quality control measures (as described above) can be masked by specifying one or more QC parameters to apply in \"sampleQCChecks\". To mask all samples failing on any quality control parameter, \"sampleQCChecks\" would be set to:\n",
    "\n",
    "```\n",
    "sampleQCChecks = ['LineWidthFail', 'CalibrationFail', 'BaselineFail', 'SolventPeakFail']\n",
    "```\n",
    "\n",
    "In this case, as only one sample is failing on 'SolventPeakFail', this sample could be masked using the following: \n",
    "\n",
    "```\n",
    "nmrData.updateMasks(filterSamples=True, sampleQCChecks=['SolventPeakFail'], filterFeatures=False)\n",
    "```\n",
    "\n",
    "Here, by additionally setting preferences with the \"sampleTypes\" and \"assayRoles\" arguments, any other samples which are not required can also be masked. In this example, we limit our dataset to study samples ('SampleType.StudySample, AssayRole.Assay') and study reference samples ('SampleType.StudyPool, AssayRole.PrecisionReference') by running 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 describe feature masking in the next section, here, we can set \"filterFeatures=False\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nmrData.updateMasks(sampleQCChecks=['SolventPeakFail'], sampleTypes=[SampleType.StudySample, SampleType.StudyPool], assayRoles=[AssayRole.Assay, AssayRole.PrecisionReference], filterFeatures=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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(nmrData, 'sample summary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mask Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For NMR datasets, typically one or more sections of the spectra known to contain unwanted or un-informative signals are removed from the data (see [Sample and Feature Masks](https://npyc-toolbox.readthedocs.io/en/latest/masks.html)).\n",
    "\n",
    "This can be done automatically using the 'updateMasks' function.\n",
    "\n",
    "The standard regions automatically masked are defined in the configuration SOP, see the 'NMRDataset Objects' table in [Built-in Configuration SOPs](https://npyc-toolbox.readthedocs.io/en/latest/configuration/builtinSOPs.html). For example, the default exclusion regions for urine are between -0.2 and 0.2 ppm (TSP) and between 4.7 and 4.9 ppm (water presaturation region).\n",
    "\n",
    "As the sample mask has been specified in the previous section, here we set \"filterSamples=False\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nmrData.updateMasks(filterSamples=False, filterFeatures=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To mask a larger region around the water presaturation signal (between 4.5 and 5.0 ppm) the data could be imported with an updated \"exclusionRegions\" parameter:\n",
    "\n",
    "```\n",
    "nmrData = nPYc.NMRDataset(rawDataPath, pulseProgram=pulseProgram, sop=sop, exclusionRegions=[[-0.2,0.2],[4.5,5.0]])\n",
    "```\n",
    "\n",
    "Also, additional regions can also be masked by using 'updateMasks' with the additional \"exclusionRegions\" parameter. For example, to additionally mask the region between 8.4 and 8.5 ppm the following would be run:\n",
    "\n",
    "```\n",
    "nmrData.updateMasks(filterSamples=False, filterFeatures=True, exclusionRegions=[(8.4, 8.5)])\n",
    "```\n",
    "\n",
    "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(nmrData, 'feature summary', withExclusions=True)"
   ]
  },
  {
   "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 the sample of unknown type with 'Sample File Name' 'DEVSET U 1D NMR raw data files/930':\n",
    "\n",
    "```\n",
    "nmrData.excludeSamples(['DEVSET U 1D NMR raw data files/930'], on='Sample File Name', message='Unknown type')\n",
    "```\n",
    "\n",
    "Or to exclude all features with 'ppm' > 8:\n",
    "\n",
    "```\n",
    "nmrData.excludeFeatures([nmrData.featureMetadata['ppm'][nmrData.featureMetadata['ppm'] > 8].values], on='ppm', message='ppm > 8')\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",
    "nmrData.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": [
    "nmrData.applyMasks()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. 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. For NMR, when smaller variables are more likely to be background noise, it may be that mean-centering the data only (\"scaling=0\") can be appropriate.\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(nmrData, 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(nmrData, 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."
   ]
  },
  {
   "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 nmrData.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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = nPYc.plotting.plotScoresInteractive(nmrData, 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 1D NMR 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(nmrData, PCAmodel, component=2)\n",
    "iplot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. 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": true
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(nmrData, '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",
    "nmrData.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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nmrData.exportDataset(saveFormat='UnifiedCSV', destinationPath='.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, to export the 'final report':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nPYc.reports.generateReport(nmrData, '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": [
    "nmrData.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",
    "nmrData.name = nmrData.name + '_PQN_normalised'\n",
    "\n",
    "# Export dataset\n",
    "nmrData.exportDataset(destinationPath='.')\n",
    "\n",
    "# Export final summary report\n",
    "nPYc.reports.generateReport(nmrData, 'final report', destinationPath='.')# Append name\n",
    "nmrData.name = nmrData.name + '_PQN_normalised'\n",
    "\n",
    "# Export dataset\n",
    "nmrData.exportDataset(destinationPath='.')\n",
    "\n",
    "# Export final summary report\n",
    "nPYc.reports.generateReport(nmrData, '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
}