ptairMS: Processing and analysis of PTR-TOF-MS data
30 June 2022
Package
ptairMS 1.4.1
1 Installation
if(!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ptairMS")2 Introduction
The ptairMS package provides a workflow to process PTR-TOF-MS raw data in the open Hierarchical Data Format 5 (HDF5; .h5 extension), and generate the peak table as an ExpressionSet object for subsequent data analysis with the many methods and packages available in R. Applications include the analysis of exhaled breath, cell culture headspace or ambient air. The package offers several features to check the raw data and tune the few processing parameters. It also enables to include new samples in a study without re-processing all the previous data, providing a convenient management for cohort studies (e.g. duration of the inclusion process, longitudinal studies, etc.).
Proton Transfer Reaction - Mass Spectrometry (PTR-MS) has emerged with excellent sensitivity and specificity for VOC analysis in a wide range of applications, including environment, food quality, and biology (Blake, Monks, and Ellis 2009). In the area of health and care, PTR-MS opens up unique opportunities for real-time analysis at the patient’s bedside.
3 Volatolomics
The characterization of volatile organic compounds (VOCs) emitted by living organisms is of major interest in medicine, food sciences, and ecology. As an example, thousands of VOCs have been identified in the exhaled breath, resulting from normal metabolism or pathological processes (Lacy Costello et al. 2014). The main advantage of breath analysis in medicine is that the sampling is non-invasive (Devillier et al. 2017). Methods based on mass spectrometry (MS) are the reference technologies for VOC analysis because of their sensitivity and large dynamic range.
4 The ptairMS processing workflow
The workflow consists of five steps:
createPtrSet: AptrSetobject is generated by taking as input the name of the directory containing the raw files (in HDF5 format), possibly grouped into subfolders according to classes of samplesdetectPeak: peak detection and quantification are performed within each file and theptrSetobject is updated with the sample metadata, the peak list for each sample, and several quality metricsalignSamples: The peak lists are aligned between samples and anExpressionSetobject is returned, containing the table of peak intensities, the sample metadata, and the feature metadata (which can be accessed with theexprs,pDataandfDatamethods from theBiobasepackage, respectively)imputing: Missing values in the table of intensities may be replaced by the integrated signal in the expected raw data regionannotateVOC: Features may be annotated, based on the Human Breathomics Database (Kuo et al. 2020)
5 Datasets (ptairData package)
Two real PTR-TOF-MS raw data sets are available in the ptairData package.
exhaledAir: Exhaled air from two healthy individual: three acquisitions per individual (with two expirations per acquisition), on distinct days [6 files]
mycobacteria: cell culture headspace: two replicates from two species and one control (culture medium) [6 files]
Note: To limit the size of the data and speed up the analysis, the raw data are truncated in the mass dimension within the range \([20.4,21.6] \cup [50,150]\) for individuals, and \([20.4,21.6] \cup [56.4,90.6]\) for mycobacteria.
6 Hands on
library(ptairMS)
library(ptairData)6.1 Graphical Interface
The whole workflow of ptairMS can be run interactively through a graphical user interface, which provides visualizations (expiration phases, peaks in the raw data, peak table, individual VOCs), quality controls (calibration, resolution, peak shape and evolution of reagent ions depending on time), and exploratory data analysis.
ptairMS::RunShinnyApp()Alternatively, the workflow can be run on the command line. We describe hereafter the main steps and parameters.
6.2 createPtrSet: Checking raw data and setting parameters
The analysis starts from a directory containing the raw data. For this example, we use the exhaledAir from the ptairData package:
dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")Before processing the data, there are two important steps:
calibrate the mass axis
determine the time periods of interest within the acquisitions (i.e., expirations or headspace analysis duration in our examples)
To perform these steps and check the quality of the raw data, we first use the createPtrSet function which, for each file:
reads the raw data
performs calibration on the
calibrationPeriodsum spectra, with themzCalibRefreference masses (see the Calibrating the mass axis section)determines the expiration time limits (or headspace duration for bacteria example) on the trace of
mzBreathTracermass or on the total sum trace if this parameter isNULL, with thefracMaxTICparameter (see the Determine the time limits for expirations or headspace duration section)builds a default
sampleMetadatatable with the file names, date and possibly the subdirectories
This function creates a unique ptrSet object for the whole study, which contains all necessary information about the files: calibration parameters, time limits, peak lists (at this step they are empty), and primary ion quantification. The corresponding raw files are available in the ptairData companion package.
Importantly, the ptrSet may be updated if new files are added or deleted from the directory, with the function updatePtrSet (see the Updating the ptrSet section).
To create a ptrSet object form the raw data directory:
exhaledPtrset <- createPtrSet(dir=dirRaw,
setName="exhaledPtrset",
mzCalibRef = c(21.022, 60.0525),fracMaxTIC = 0.7,
saveDir = NULL )## ind1-1.h5 check## ind1-2.h5 check## ind1-3.h5 check## ind2-1.h5 check## ind2-2.h5 check## ind2-3.h5 checkexhaledPtrset## ptrSet object : exhaledPtrset
## directory: /home/biocbuild/bbs-3.15-bioc/R/library/ptairData/extdata/exhaledAir
## 6 files contained in the directory
## 6 files checked
## 0 files processed by detectPeakTo get the list of the file names in your directory (at the last time the ptrSet was created or updated), use getFileNames:
getFileNames(exhaledPtrset)## [1] "ind1-1.h5" "ind1-2.h5" "ind1-3.h5" "ind2-1.h5" "ind2-2.h5" "ind2-3.h5"To check the quality and view useful information about your files, you can use the plot method on the ptrSet, that provides four plots:
calibError: the calibration error boxplot in ppm\(\frac{m}{\Delta m}\): an estimation of the resolution for the
mzCalibRefpeaksthe average normalized shape of the
mzCalibRefpeaksThe primary ion isotope intensity in cps in function of acquisition date
plot(exhaledPtrset)
We will now quickly explain each step of this function, and show how to choose and check the quality of the mzCalibRef, mzBreathTracer and fracMaxTIC parameters.
6.2.1 Calibrating the mass axis
To convert Time Of Flight (TOF) axis to mass-to-charge ratio (mz), the following formula is used (Müller et al. 2013):
- \(mz = \Big (\frac{TOF-b}{a} \Big)^2\)
To estimate the parameters (a,b), reference peaks with known masses (and without overlapping peaks) are needed. For optimal results, the calibration peaks have to be well distributed over the entire m/z axis. The masses of the calibration peaks are set with the mzCalibRef parameter, in an numeric vector. For exhaled air, we propose by default 6 peaks (two of them are suggested by the IONICON manufacturer):
calib_table<-read.csv(system.file("extdata", "reference_tables/calib_table.tsv", package = "ptairMS"),sep="\t")
knitr::kable(calib_table)| Formula | mz | Compound |
|---|---|---|
| H3 (18O)+ | 21.02204 | Isotope of primary ion (H3O+) |
| N2H+ | 29.01342 | Diazote |
| C3H5+ | 41.03858 | Allyl cation |
| C3H7O+ | 60.05250 | Isotope of Acetone |
| C6H6I+ | 203.94300 | Iodobenzene (1,3-diiodobenzene fragment) |
| C6H5I2+ | 330.84950 | 1,3-diiodobenzene (instrument internal standard) |
Other reference masses (at least two) can be provided by the user.
Since there may be a drift of mass calibration of the instrument during the acquisition, e.g. due to low change of temperature, periodic calibration is performed every calibrationPeriods seconds (default 60): correction of the drifts is then performed by linear interpolation, by using the first calibration as the reference.
If a mass contains one or several outlier values on the summary plot (i.e. an error > 20 ppm) for specific files, you can check those files with the function plotCalib that plots the Average total ion spectrum around all the mzCalibRef (with the exact masses as red vertical lines), after calibration drift correction.
plotCalib(exhaledPtrset,fileNames=getFileNames(exhaledPtrset)[1])
You can then choose to keep or delete these masses from mzCalibRef. To change the calibration masses, use the calibration function on the ptrSet object.
In our dataset example, since we have truncated the mass axis, only two masses from those listed above are present. To calibrate with three masses, we therefore add the peak 75.04406 ([C3H6O2+H]+, Hydroxyacetone).
exhaledPtrset <- calibration(exhaledPtrset, mzCalibRef = c(21.022, 60.0525,75.04406))
plot(exhaledPtrset,type="calibError")
6.2.2 Determine the time limits of expirations or headspace duration
To see and change the expirations or headspace duration, use the changeTimeLimits function that opens the Shiny application for interactive visualization and selection, or use the timeLimits function.
The limits are determined on the total trace if mzBreathTracer is NULL, or on the trace around mzBreathTracer mass. They correspond to the part of the trace where the intensity is higher than fracMaxTIC * max(TIC), after baseline removal if baseline is set to TRUE.
Note: To analyze the entire spectrum (e.g. in ambient air studies), set fracMaxTIC = 0.
exhaledPtrset <- changeTimeLimits(exhaledPtrset)Example of expiration detection at fracMaxTIC = 0.5 on a single file:
samplePath <-getFileNames(exhaledPtrset,fullNames = TRUE)[1]
sampleRaw <- readRaw(samplePath, calib = FALSE)
expirationLimit <- timeLimits(sampleRaw,fracMaxTIC = 0.5,plotDel = TRUE, mzBreathTracer = 60.05)
expirationLimit <- timeLimits(sampleRaw,fracMaxTIC = 0.9,plotDel = TRUE,mzBreathTracer = NULL)
Note: You can also see the expiration limits on all or several files from the ptrSet with the plotTIC function. By default, when fileNames is NULL, all TICs files are plotted. A pdf file may be generated by setting the pdfFile argument to the absolute file path ending with the .pdf extension. Finally, you can remove the baseline by setting baselineRm = TRUE, and add the time limits to the plot by setting showLimits = TRUE. Coloring the TICs according to a column from the sampleMetadata is also possible by indicating the column name as the colorBy parameter.
plotTIC(object = exhaledPtrset,baselineRm = TRUE,type = "ggplot")
6.2.3 Managing sample metadata
The createPtrSet function automatically generates a default sampleMetadata data.frame. It contains the file names as row names, a column named ‘subfolder’ when the files are organized into subfolders in the parent directory, and a column date with the date and hour of the acquisition. To get this data frame, use the getSampleMetadata method. Remember that the row names of the sampleMetadata must always correspond to all the file names from the directory.
getSampleMetadata(exhaledPtrset)## subfolder date
## ind1-1.h5 ind1 21/02/2019 11:52:49
## ind1-2.h5 ind1 27/02/2019 11:41:24
## ind1-3.h5 ind1 28/02/2019 11:28:29
## ind2-1.h5 ind2 07/02/2019 10:29:33
## ind2-2.h5 ind2 18/02/2019 11:06:52
## ind2-3.h5 ind2 22/02/2019 09:52:20You can at any moment obtain this default sample metadata with the function resetSampleMetadata(exhaledPtrset).
To modify the sample metadata, there exists two different ways:
setSampleMetadatamethod: modify the data frame in your R session and set it back into theptrSetobject.
sampleMD <- getSampleMetadata(exhaledPtrset)
colnames(sampleMD)[1] <- "individual"
exhaledPtrset <- setSampleMetadata(exhaledPtrset,sampleMD)
getSampleMetadata(exhaledPtrset)## individual date
## ind1-1.h5 ind1 21/02/2019 11:52:49
## ind1-2.h5 ind1 27/02/2019 11:41:24
## ind1-3.h5 ind1 28/02/2019 11:28:29
## ind2-1.h5 ind2 07/02/2019 10:29:33
## ind2-2.h5 ind2 18/02/2019 11:06:52
## ind2-3.h5 ind2 22/02/2019 09:52:20exportSampleMetadaandimportSampleMetadatamethods:exportSampleMetadasaves the data.frame in the.tsvformat in the directory of your choice. The parametersaveFilemust always have the.tsvextension. You can open and modify it on any spreadsheet editor, but row names must never be modified: if one of the files happens to be missing from the row names during the subsequent import, an error will be thrown. Once the.tsvfile has been modified, it can be imported back into theptrSetobject, with the functionimportSampleMetadata.
exportSampleMetada(exhaledPtrset, saveFile = file.path(DirBacteria,"sampleMetadata.tsv"))
exhaledPtrset <- importSampleMetadata(exhaledPtrset, file = file.path(DirBacteria,"sampleMetadata.tsv"))6.2.4 Saving
When calling the createPtrSet function, you may use the saveDir argument to save the ptrSet object in the directory of your choice, with setName parameter as name (the .RData extension will automatically be added at the end of the file name). Subsequent import of the saved ptrSet object relies on the classical load function.
6.2.5 Plot raw data
There exist two functions for plotting the raw data:
plotRaw: Displays for a file the image of the matrix of intensities, the Extracted Ion Chromatogram (EIC), the Extracted Ion Spectrum (EIS), for the selected m/z and time ranges, and all the putative annotations available in thevocDBwithin this mz range (no peak detection has yet been done as this step)
plotRaw(exhaledPtrset, mzRange = 59 , fileNames = getFileNames(exhaledPtrset)[1],showVocDB = TRUE)
## ion_mass ion_formula name_iupac
## 40 59.08552 [C4H10+H]+ butane|2-methylpropane
## 39 59.06037 [C2H6N2+H]+ dimethyldiazene
## 38 59.04914 [C3H6O+H]+ prop-2-en-1-ol|propanal|propan-2-one
## 37 58.94261 [Ni+H]+ nickelplotFeatures: for all selected files, the average spectrum for all time periods is plotted, with the background superimposed as a dashed line. As with theplotTICfunction, you can choose the type of display (plotlyorggplot), and the possibility to color spectra individually by indicating a column name of the sample Metadata
plotFeatures(exhaledPtrset,mz=59.049,type="ggplot",colorBy = "individual")
6.3 updatePtrSet: Updating the ptrSet
If you delete or add files to the directory after the ptrSet object has been created, the updatePtrSet function must be run.
exhaledPtrset <- updatePtrSet(exhaledPtrset)## nothing to update6.4 detectPeak: Peak detection and quantification
Now that we have checked that the calibration is efficient, and that the expiration or headspace time limits are correct, the peaks can be detected and quantified in each file with the detectPeak function, which works on the ptrSet object (and the corresponding raw files).
For each file in the directory, this function:
performs the calibration of the mass axis every
calibrationPeriodseconddetects peaks on the total ion spectrum with an untargeted peak picking algorithm
estimates the temporal evolution of each detected peak with a 2-dimensional model
quantifies each peak in both exhaled air and background, in cps, ncps (normalized by the primary ion, H3(O18)+ * 488) and ppb [(Hansel et al. 1995)] units, and stores the corresponding information in the
peakListslot of the ptrSet object as a list of ExpressionSet objects (see below)performs two unilateral t-test comparing the mean intensities between background and expiration/headspace phases, and stores the p-value in the
peakListAlignedslot
The peakList is then written in the ptrSet object as a list of ExpressionSet, each containing all the peaks detected in one file (e.g., sample; see below for the details of the ExpressionSet content).
exhaledPtrset <- detectPeak(exhaledPtrset)## ind1-1.h5 : 157 peaks detected
## ind1-2.h5 : 168 peaks detected
## ind1-3.h5 : 141 peaks detected
## ind2-1.h5 : 166 peaks detected
## ind2-2.h5 : 177 peaks detected
## ind2-3.h5 : 163 peaks detectedThe peak detection may be restricted to specific nominal masses with the mzNominal argument: for example exhaledPtrset <- detectPeak(exhaledPtrset , mzNominal = c(5,60)).
The peak detection step may take a few minutes if there are many files and a large m/z range (1 to 2 minutes for files for files with an average acquisition time of 3 minutes) . Parallel computing is available by setting parallelize = TRUE and by giving the number of available cores of your computer in nbCores. To find the number of CPU cores available in your computer, use parallel::detectCores().
To see the resulting peak lists, use the getPeakList method. It returns a list of ExpressionSet, where:
assay Data: the matrix of peak intensities, with m/z peak centre as row names, and the quantification in cps at each time point
feature Data: the data frame with m/z peak centre as row names, and the following columns:
- quanti_unit: the mean of the quantifications for all expiration/headspace phases
- background_unit: the mean of the quantifications for all background phases
- diffAbs_unit: the mean of the quantifications for all expiration/headspace phases after baseline subtraction
- pValLess/ pValGreater: the p-value of the unilateral t-test, which tests that quantification (in cps) of expiration points are less/greater than the intensity of the background.
- lower/upper: integration boundaries
- parameter peak: estimated peak parameters (width et left and rigth and height)
peakList<-getPeakList(exhaledPtrset)
peakList1<-peakList$`ind1-1.h5`
X<-Biobase::exprs(peakList1)
Y<-Biobase::fData(peakList1)
mz<-Y[,"Mz"]
plot(X[which.min(abs(mz-59.0498)),],ylab="cps",xlab="time",main=paste("Temporal evolution of acetone "))
head(Y)## Mz lowerMz upperMz overlap parameterPeak.delta1
## 21.022 21.02203 21.01170 21.03136 0 0.002024031
## 51.0443 51.04429 51.02268 51.06376 0 0.005305088
## 53.0032 53.00323 52.97378 53.02971 1 0.006625404
## 53.0389 53.03889 53.01388 53.06145 1 0.004634475
## 54.0077 54.00772 53.97772 54.03470 1 0.006750965
## 54.0362 54.03616 54.00614 54.06315 1 0.006754520
## parameterPeak.delta2 parameterPeak.height quanti_cps background_cps
## 21.022 0.002627754 819.62724 3955.9945 4619.7859
## 51.0443 0.004436511 443.84840 5538.3828 122.4474
## 53.0032 0.006625404 101.43480 691.3785 1404.7156
## 53.0389 0.006629563 89.39703 756.2871 709.5341
## 54.0077 0.006750965 20.57361 196.0334 246.1388
## 54.0362 0.006754520 17.89411 124.4829 260.2830
## diffAbs_cps pValLess pValGreater quanti_ncps background_ncps
## 21.022 600.5037 2.553555e-13 1.000000e+00 1.886460e-03 0.0022029958
## 51.0443 5366.5027 1.000000e+00 4.145078e-23 2.641039e-03 0.0000583904
## 53.0032 813.5105 2.923887e-11 1.000000e+00 3.296914e-04 0.0006698541
## 53.0389 218.7715 1.000000e+00 1.378881e-01 3.606438e-04 0.0003383491
## 54.0077 16.8247 7.036460e-04 1.000000e+00 9.348067e-05 0.0001173740
## 54.0362 146.0940 1.322352e-09 1.000000e+00 5.936102e-05 0.0001241188
## diffAbs_ncps quanti_ppb background_ppb diffAbs_ppb
## 21.022 2.863568e-04 45.09934453 52.34038369 6.803473751
## 51.0443 2.559076e-03 2.80151579 0.06155455 2.697751545
## 53.0032 3.879314e-04 0.33538163 0.67719280 0.392181478
## 53.0389 1.043236e-04 0.36686824 0.34205597 0.105466553
## 54.0077 8.023044e-06 0.09318326 0.11627559 0.007947961
## 54.0362 6.966652e-05 0.05917216 0.12295727 0.0690145576.5 Updating the ptrSet peak lists with detectPeak
As described previously, an important feature of the ptairMS package is the possibility to update the ptrSet object linked to a directory, as new data files are added and included in this directory. In such a case, we have seen that the updatePtrSet function must be used to reset the sample metadata or append it. Then, the detectPeak function must be used on the updated ptrSet to compute the additional peak lists.
exhaledPtrset<-updatePtrSet(exhaledPtrset)
exhaledPtrset<-setSampleMetadata(exhaledPtrset,resetSampleMetadata(exhaledPtrset))
exhaledPtrset<-detectPeak(exhaledPtrset)6.6 alignSamples: Aligning features between samples
The alignment between samples (i.e. the matching of variables between the peak lists within the ptrSet object) is performed by using a kernel gaussian density (Delabrière et al. 2017).
The alignSamples returns an ExpressionSet object, with the table of peak intensities which has just been built, the sample meta data (borrowed from the input ptrSet) and the variable meta data which contains peak intensities in the background.
It is possible to apply two filters:
On the background: only variables with a significant higher intensity in the expiration phases compared to the background phases (at the
pValGreaterThresp-value threshold) forfracExp% of the samples are keptOn sample reproducibility: only variables which are detected in more than
fracGrouppercent of the samples (or percent of at least onegroup) are kept.
If you do not want to apply those filters, set fracGroup to 0 and pValGreaterThres to 1.
exhaledEset <- alignSamples(exhaledPtrset, group="individual", fracGroup = 1, fracExp=1/6)## 113 peaks alignedThe three tables from the ExpressionSet can be accessed with the classical exprs, pData, and fData accessors:
knitr::kable(head(Biobase::exprs(exhaledEset)))| ind1-1.h5 | ind1-2.h5 | ind1-3.h5 | ind2-1.h5 | ind2-2.h5 | ind2-3.h5 | |
|---|---|---|---|---|---|---|
| 51.0442 | 2.8015158 | 2.7043680 | 2.3562760 | 8.1931241 | 6.9344951 | 10.9064801 |
| 52.0473 | NA | NA | NA | 0.1242310 | 0.1099690 | 0.1572898 |
| 53.0028 | 0.3353816 | 0.4038502 | 0.3273306 | 0.3843748 | 0.6015839 | 0.4769054 |
| 53.0137 | NA | NA | NA | 0.1215173 | 0.1738785 | 0.1656050 |
| 53.0388 | 0.3668682 | 0.3506322 | 0.3529729 | 0.5996653 | 0.7446162 | 0.7024514 |
| 54.0086 | 0.0931833 | 0.0933895 | NA | 0.1436999 | 0.2591290 | 0.1908737 |
knitr::kable(head(Biobase::pData(exhaledEset)))| individual | date | |
|---|---|---|
| ind1-1.h5 | ind1 | 21/02/2019 11:52:49 |
| ind1-2.h5 | ind1 | 27/02/2019 11:41:24 |
| ind1-3.h5 | ind1 | 28/02/2019 11:28:29 |
| ind2-1.h5 | ind2 | 07/02/2019 10:29:33 |
| ind2-2.h5 | ind2 | 18/02/2019 11:06:52 |
| ind2-3.h5 | ind2 | 22/02/2019 09:52:20 |
knitr::kable(head(Biobase::fData(exhaledEset)))| ion_mass | Background - ind1-1.h5 | Background - ind1-2.h5 | Background - ind1-3.h5 | Background - ind2-1.h5 | Background - ind2-2.h5 | Background - ind2-3.h5 | |
|---|---|---|---|---|---|---|---|
| 51.0442 | 51.0442 | 0.0615546 | 0.0368512 | 0.0364317 | 0.0689280 | 0.0230780 | 0.1019756 |
| 52.0473 | 52.0473 | NA | NA | NA | 0.0089867 | 0.0056339 | 0.0096605 |
| 53.0028 | 53.0028 | 0.6771928 | 0.4488869 | 0.1985743 | 0.9340833 | 0.4764598 | 0.6112486 |
| 53.0137 | 53.0137 | NA | NA | NA | 0.1573156 | 0.0474827 | 0.1436644 |
| 53.0388 | 53.0388 | 0.3420560 | 0.3573706 | 0.2133881 | 0.1857090 | 0.3224024 | 0.3629669 |
| 54.0086 | 54.0086 | 0.1162756 | 0.0794680 | NA | 0.2342243 | 0.0949972 | 0.1240503 |
Note: The view method from the ropls package may be used to print these tables:
6.7 imputing: Imputation of missing values
To impute missing values, ptairMS returns back to the raw data, and performs the quantification with the same method as detectPeak but this time without any limit of detection.
exhaledEset <- ptairMS::impute(exhaledEset, exhaledPtrset)## ind1-1.h5 done## ind1-2.h5 done## ind1-3.h5 done## ind2-1.h5 done## ind2-2.h5 done## ind2-3.h5 done6.8 annotateVOC: Annotation
ptairMS provides putative annotations by matching the measured ion masses to the Human Breathomics Database (Kuo et al. 2020). Applied to an ExpressionSet, it appends the feature metadata (fData) with new columns containing chemical information (formulas, IUPAC name, InChI, etc.), as well as the isotopes for nuclides C13, O17 and O18.
annotateVOC(59.049)## vocDB_ion_mass vocDB_ion_formula vocDB_name_iupac
## 59.049 59.04914 [C3H6O+H]+ prop-2-en-1-ol, propan-2-one, propanalexhaledEset<-annotateVOC(exhaledEset)
knitr::kable(head(Biobase::fData(exhaledEset)))| ion_mass | Background - ind1-1.h5 | Background - ind1-2.h5 | Background - ind1-3.h5 | Background - ind2-1.h5 | Background - ind2-2.h5 | Background - ind2-3.h5 | vocDB_ion_mass | vocDB_ion_formula | vocDB_name_iupac | isotope | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 51.0442 | 51.0442 | 0.0615545543 | 0.0368512260 | 0.0364317310 | 0.0689280372 | 0.0230780111 | 0.1019755646 | 53.0388/52.0473 | |||
| 52.0473 | 52.0473 | NA | NA | NA | 0.0089867492 | 0.0056338837 | 0.0096605388 | NA | |||
| 53.0028 | 53.0028 | 0.6771928040 | 0.4488868914 | 0.1985743123 | 0.9340832591 | 0.4764597682 | 0.6112485766 | 54.0086 | |||
| 53.0137 | 53.0137 | NA | NA | NA | 0.1573156277 | 0.0474827179 | 0.1436643932 | NA | |||
| 53.0388 | 53.0388 | 0.3420559656 | 0.3573705533 | 0.2133881456 | 0.1857090164 | 0.3224023506 | 0.3629669359 | 53.03857 | [C4H4+H]+ | but-1-en-3-yne | NA |
| 54.0086 | 54.0086 | 0.1162755866 | 0.0794679724 | NA | 0.2342242578 | 0.0949972399 | 0.1240502983 | NA |
6.9 writeEset: Export data and metadata to 3 tabular files
Finally, the ExpressionSet can be exported to 3 tabulated files ‘dataMatrix.tsv’, sampleMetadata.tsv’, and ‘variableMetadata.tsv’:
writeEset(exhaledEset, dirC = file.path(getwd(), "processed_dataset"))6.10 Statistical analysis
The ExpressionSet object is now ready for subsequent data analysis (e.g. data mining, classification or feature selection) with the many R packages.
As an example, we describe here how to perform Exploratory Data Analysis with the ropls package (Thevenot 2015).
As a preliminary step, log transformation is often used in Mass Spectrometry to stabilize the variance:
Biobase::exprs(exhaledEset) <- log2(Biobase::exprs(exhaledEset))Then, the data and metadata may be printed and plotted with the view method:
ropls::view(Biobase::exprs(exhaledEset),printL=FALSE)
To avoid redundancy, the isotopes may be discarded:
isotopes<-Biobase::fData(exhaledEset)[,"isotope"]
isotopes<-isotopes[!is.na(isotopes)]
exhaledEset <- exhaledEset[!(Biobase::fData(exhaledEset)[, "ion_mass"] %in% isotopes), ]Principal Component Analysis may then be performed (due to the limited number of samples, cross-validation is decreased to 5 instead of 7):
exhaledPca<-ropls::opls(exhaledEset,crossvalI=5,info.txtC="none",fig.pdfC="none")
ropls::plot(exhaledPca, parAsColFcVn=Biobase::pData(exhaledEset)[, "individual"],typeVc="x-score")
Subsequent supervised analysis, e.g. (Orthogonal) Partial Least Squares, or feature selection may be performed with the ropls and biosigner packages, respectively.
To get the most important variable that contribute to the dimension that discriminate the two individual (here the first dimension), we look at the loadings:
load1<-ropls::getLoadingMN(exhaledPca)[,1]
barplot(sort(abs(load1),decreasing = TRUE))
knitr::kable(Biobase::fData(exhaledEset)[names(sort(abs(load1),decreasing = TRUE)[1:10]),c("vocDB_ion_formula","vocDB_name_iupac")])| vocDB_ion_formula | vocDB_name_iupac | |
|---|---|---|
| 110.9703 | ||
| 53.0388 | [C4H4+H]+ | but-1-en-3-yne |
| 67.0543 | [C5H6+H]+ | (E)-pent-3-en-1-yne, 2-methylbut-1-en-3-yne, cyclopenta-1,3-diene, pent-3-en-1-yne |
| 109.0656 | [C7H8O+H]+ | 2-prop-2-enylfuran, 3-methylphenol, 4-methylphenol, phenylmethanol |
| 117.0556 | ||
| 52.0473 | ||
| 69.0704 | [C5H8+H]+ | (3E)-penta-1,3-diene, (3Z)-penta-1,3-diene, 2-methylbuta-1,3-diene, 3-methylbuta-1,2-diene, pent-2-yne, penta-1,2-diene, penta-1,4-diene |
| 51.0442 | ||
| 135.1133 | [C10H14+H]+ | 1,2,3,4-tetramethylbenzene, 1,2,3,5-tetramethylbenzene, 1,2,4,5-tetramethylbenzene, 1,3-diethylbenzene, 1,4-diethylbenzene, 1-ethyl-2,3-dimethylbenzene, 1-ethyl-2,4-dimethylbenzene, 1-ethyl-3,5-dimethylbenzene, 1-methyl-2-propan-2-ylbenzene, 1-methyl-2-propylbenzene, 1-methyl-3-propan-2-ylbenzene, 1-methyl-4-propan-2-ylbenzene, 1-methyl-4-propylbenzene, 2-ethyl-1,4-dimethylbenzene, 2-methyl-5-prop-1-en-2-ylcyclohexa-1,3-diene, 4-ethyl-1,2-dimethylbenzene, butylbenzene |
| 129.127 | [C8H16O+H]+ | 2-ethylhexanal, 3-methylheptan-2-one, 6-methylheptan-2-one, 6-methylheptan-3-one, octan-2-one, octan-3-one, octanal |
We then could plot the raw data around this masses to check the robustness of these potential marker:
plotFeatures(exhaledPtrset,mz = 53.0387,typePlot = "ggplot",colorBy = "individual")
plotFeatures(exhaledPtrset,mz = 67.0539,typePlot = "ggplot",colorBy = "individual")
6.11 Mycobacteria dataset
The package can also be used to process data from headspace analysis (e.g., cell culture headspace). The acquisition phase is detected on the total ion trace, in a similar way as for exhaled air. The mycobacteria dataset contains two replicates from two bacterial species and one control (culture medium). The 6 corresponding raw files are available in the companion ptairData package.
dir <- system.file("extdata/mycobacteria", package = "ptairData")
mycobacteriaSet <- createPtrSet(dir = dir, setName = "test",
mzCalibRef = c(21.022,59.049))## Control1.h5 check## Control2.h5 check## Specie-a1.h5 check## Specie-a2.h5 check## specie-b1.h5 check## specie-b2.h5 checkmycobacteriaSet <- detectPeak(mycobacteriaSet, smoothPenalty = 0)## Control1.h5 : 36 peaks detected
## Control2.h5 : 38 peaks detected
## Specie-a1.h5 : 38 peaks detected
## Specie-a2.h5 : 39 peaks detected
## specie-b1.h5 : 51 peaks detected
## specie-b2.h5 : 41 peaks detectedplotTIC(mycobacteriaSet,type="ggplot",showLimits = TRUE,file="Specie-a2.h5")
eSet <- alignSamples(mycobacteriaSet)## 30 peaks alignedeSet<-impute(eSet,ptrSet = mycobacteriaSet)## specie-b1.h5 done## specie-b2.h5 doneX<-Biobase::exprs(eSet)
pca<-ropls::opls(log2(t(X)),predI =2,crossvalI=5,info.txtC = "none",
fig.pdfC = "none")
ropls::view(log2(X),printL=FALSE)
plot(pca,type="x-score")
6.12 Processing a single raw file
To read and access the entire raw data of a single file, use the readRaw function. It opens the data from an absolute file path with the rhdf5 library, and returns a ptrRaw object containing the raw data matrix, the time axis, the m/z axis obtained after calibration if calibTIS = TRUE, and additional information contained in the h5 file (transmission curve, drift temperature, …).
As an example, we get the absolute file path for the first raw file:
dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
samplePath<-list.files(dirRaw,recursive = TRUE,full.names = TRUE,pattern = ".h5$")[1]To read the file:
sampleRaw <- readRaw(samplePath, calib = FALSE)
sampleRaw## /home/biocbuild/bbs-3.15-bioc/R/library/ptairData/extdata/exhaledAir/ind1/ind1-1.h5
## mz range: 20.4 - 150.7
## time range: 0 - 49
## Calibration error in ppm:
## No calibration performsAt this stage, the same methods as with the ptrSet object can be used: calibration, timeLimits, plotCalib, plotTIC, plotRaw, and detectPeak.
7 Acknowledgements
This research was conducted by Camille Roquencourt, with the contributions from Paul Zheng, Pierrick Roger-Mele and Etienne A. Thevenot (Data Sciences for Deep Phenotyping and Precision Medicine team at CEA), in collaboration with Stanislas Grassin-Delyle, Helene Salvator, Emmanuel Naline, Philippe Devillier and Louis-Jean Couderc (Exhalomics platform at the Foch Hospital) within the SoftwAiR project funded by the Agence Nationale de la Recherche (ANR-18-CE45-0017 grant).
8 Session Info
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ptairData_1.4.0 ptairMS_1.4.1 BiocStyle_2.24.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 ggsignif_0.6.3
## [3] ellipsis_0.3.2 htmlTable_2.4.0
## [5] qqman_0.1.8 XVector_0.36.0
## [7] GenomicRanges_1.48.0 base64enc_0.1-3
## [9] rstudioapi_0.13 ggpubr_0.4.0
## [11] farver_2.1.0 MultiDataSet_1.24.0
## [13] MultiAssayExperiment_1.22.0 bit64_4.0.5
## [15] fansi_1.0.3 codetools_0.2-18
## [17] splines_4.2.0 doParallel_1.0.17
## [19] knitr_1.39 Formula_1.2-4
## [21] jsonlite_1.8.0 broom_0.8.0
## [23] cluster_2.1.3 png_0.1-7
## [25] shinyscreenshot_0.2.0 shiny_1.7.1
## [27] BiocManager_1.30.18 compiler_4.2.0
## [29] backports_1.4.1 assertthat_0.2.1
## [31] Matrix_1.4-1 fastmap_1.1.0
## [33] limma_3.52.2 cli_3.3.0
## [35] later_1.3.0 htmltools_0.5.2
## [37] tools_4.2.0 GenomeInfoDbData_1.2.8
## [39] gtable_0.3.0 glue_1.6.2
## [41] dplyr_1.0.9 Rcpp_1.0.8.3
## [43] carData_3.0-5 Biobase_2.56.0
## [45] jquerylib_0.1.4 vctrs_0.4.1
## [47] rhdf5filters_1.8.0 iterators_1.0.14
## [49] xfun_0.31 stringr_1.4.0
## [51] enviPat_2.4 mime_0.12
## [53] lifecycle_1.0.1 rstatix_0.7.0
## [55] zlibbioc_1.42.0 MASS_7.3-57
## [57] scales_1.2.0 MatrixGenerics_1.8.1
## [59] promises_1.2.0.1 SummarizedExperiment_1.26.1
## [61] parallel_4.2.0 rhdf5_2.40.0
## [63] RColorBrewer_1.1-3 ropls_1.28.2
## [65] yaml_2.3.5 gridExtra_2.3
## [67] ggplot2_3.3.6 sass_0.4.1
## [69] rpart_4.1.16 calibrate_1.7.7
## [71] latticeExtra_0.6-29 stringi_1.7.6
## [73] highr_0.9 S4Vectors_0.34.0
## [75] foreach_1.5.2 checkmate_2.1.0
## [77] BiocGenerics_0.42.0 chron_2.3-57
## [79] GenomeInfoDb_1.32.2 matrixStats_0.62.0
## [81] rlang_1.0.3 pkgconfig_2.0.3
## [83] bitops_1.0-7 evaluate_0.15
## [85] lattice_0.20-45 purrr_0.3.4
## [87] Rhdf5lib_1.18.2 htmlwidgets_1.5.4
## [89] labeling_0.4.2 cowplot_1.1.1
## [91] bit_4.0.4 tidyselect_1.1.2
## [93] magrittr_2.0.3 bookdown_0.27
## [95] R6_2.5.1 IRanges_2.30.0
## [97] magick_2.7.3 generics_0.1.2
## [99] Hmisc_4.7-0 DelayedArray_0.22.0
## [101] DBI_1.1.3 pillar_1.7.0
## [103] foreign_0.8-82 survival_3.3-1
## [105] abind_1.4-5 RCurl_1.98-1.7
## [107] nnet_7.3-17 tibble_3.1.7
## [109] crayon_1.5.1 car_3.1-0
## [111] utf8_1.2.2 rmarkdown_2.14
## [113] jpeg_0.1-9 grid_4.2.0
## [115] minpack.lm_1.2-2 data.table_1.14.2
## [117] digest_0.6.29 xtable_1.8-4
## [119] tidyr_1.2.0 httpuv_1.6.5
## [121] signal_0.7-7 stats4_4.2.0
## [123] munsell_0.5.0 bslib_0.3.1Bibliography
Blake, Robert S., Paul S. Monks, and Andrew M. Ellis. 2009. “Proton-Transfer Reaction Mass Spectrometry.” Chemical Reviews 109 (3): 861–96. https://doi.org/10.1021/cr800364q.
Delabrière, Alexis, Ulli M Hohenester, Benoit Colsch, Christophe Junot, François Fenaille, and Etienne A Thévenot. 2017. “proFIA: A Data Preprocessing Workflow for Flow Injection Analysis Coupled to High-Resolution Mass Spectrometry.” Edited by Oliver Stegle. Bioinformatics 33 (23): 3767–75. https://doi.org/10.1093/bioinformatics/btx458.
Devillier, Philippe, Helene Salvator, Emmanuel Naline, Louis-Jean Couderc, and Stanislas Grassin-Delyle. 2017. “Metabolomics in the Diagnosis and Pharmacotherapy of Lung Diseases.” Current Pharmaceutical Design 23 (14). https://doi.org/10.2174/1381612823666170130155627.
Hansel, A., A. Jordan, R. Holzinger, P. Prazeller, W. Vogel, and W. Lindinger. 1995. “Proton Transfer Reaction Mass Spectrometry: On-Line Trace Gas Analysis at the Ppb Level.” International Journal of Mass Spectrometry and Ion Processes 149-150 (November): 609–19. https://doi.org/10.1016/0168-1176(95)04294-U.
Kuo, Tien-Chueh, Cheng-En Tan, San-Yuan Wang, Olivia A Lin, Bo-Han Su, Ming-Tsung Hsu, Jessica Lin, et al. 2020. “Human Breathomics Database” 2020: 8.
Lacy Costello, B de, A Amann, H Al-Kateb, C Flynn, W Filipiak, T Khalid, D Osborne, and N M Ratcliffe. 2014. “A Review of the Volatiles from the Healthy Human Body.” Journal of Breath Research 8 (1): 014001. https://doi.org/10.1088/1752-7155/8/1/014001.
Müller, Markus, Tomáš Mikoviny, Werner Jud, Barbara D’Anna, and Armin Wisthaler. 2013. “A New Software Tool for the Analysis of High Resolution PTR-TOF Mass Spectra.” Chemometrics and Intelligent Laboratory Systems 127 (August): 158–65. https://doi.org/10.1016/j.chemolab.2013.06.011.
Thevenot. 2015. “Analysis of the Human Adult Urinary Metabolome Variations with Age, Body Mass Index, and Gender by Implementing a Comprehensive Workflow for Univariate and Opls Statistical Analyses.”