Chapter 6 Raw data pretreatment
Raw data from the instruments such as LC-MS or GC-MS were hard to be analyzed. To make it clear, the structure of those data could be summarized as:
Indexed scan with time-stamp
Each scan contains a full scan mass spectra
Common formats for open source mass spectrum data are mzxml, mzml or CDF. However, MassComp might shrink the data size(R. Yang, Chen, and Ochoa 2019).
mzML2ISA & nmrML2ISA could generate enriched ISA-Tab metadata files from metabolomics XML data (Larralde et al. 2017).
6.1 Data visualization
6.2 Peak extraction
GC/LC-MS data are usually be shown as a matrix with column standing for retention times and row standing for masses after bin them into small cell.
Conversation from the mass-retention time matrix into a vector with selected MS peaks at certain retention time is the basic idea of Peak extraction. You could EIC for each mass to charge ratio and use the change of trace slope to determine whether there is a peak or not. Then we could make integration for this peak and get peak area and retention time.
<- c(10,10,10,10,10,14,19,25,30,33,26,21,16,12,11,10,9,10,11,10) intensity <- c(1:20) time plot(intensity~time, type = 'o', main = 'EIC')
However, due to the accuracy of instrument, the detected mass to charge ratio would have some shift and EIC would fail if different scan get the intensity from different mass to charge ratio.
matchedfilter algorithm (Smith et al. 2006), they solve this issue by bin the data in m/z dimension. The adjacent chromatographic slices could be combined to find a clean signal fitting fixed second-derivative Gaussian with full width at half-maximum (fwhm) of 30s to find peaks with about 1.5-4 times the signal peak width. The the integration is performed on the fitted area.
Centwave algorithm (Tautenhahn, Böttcher, and Neumann 2008) based on detection of regions of interest(ROI) and the following Continuous Wavelet Transform (CWT) is preferred for high-resolution mass spectrum. ROI means a region with stable mass for a certain time. When we find the ROIs, the peak shape is evaluated and ROI could be extended if needed. This algorithm use
prefilter to accelerate the processing speed.
prefilter with 3 and 100 means the ROI should contain 3 scan with intensity above 100. Centwave use a peak width range which should be checked on pool QC. Another important parameter is
ppm. It is the maximum allowed deviation between scans when locating regions of interest (ROIs), which is different from vendor number and you need to extend them larger than the company claimed. For
profparam, it’s used for fill peaks or align peaks instead of peak picking.
snthr is the cutoff of signal to noise ratio.
Deep learning frame for LC-MS feature detection on 2D pseudo color image could improve the peak picking process (F. Zhao, Huang, and Zhang 2021).
Various data acquisition workflow could be checked here(Fenaille et al. 2017b).
decoMS2 An Untargeted Metabolomic Workflow to Improve Structural Characterization of Metabolites(Nikolskiy et al. 2013). It requires two different collision energies, low (usually 0V) and high, in each precursor range to solve the mathematical equations.
Data-Independent Targeted Metabolomics Method could connect MS1 and MRM (Y. Chen et al. 2017)
The coverage of DDA could be enhanced by a feature classification strategy (Hu, Cai, and Huan 2019).
MetaboDIA workflow build customized MS/MS spectral libraries using a user’s own data dependent acquisition (DDA) data and to perform MS/MS-based quantification with DIA data, thus complementing conventional MS1-based quantification (G. Chen et al. 2017)
Skyline is a freely-available and open source Windows client application for building Selected Reaction Monitoring (SRM) / Multiple Reaction Monitoring (MRM), Parallel Reaction Monitoring (PRM - Targeted MS/MS), Data Independent Acquisition (DIA/SWATH) and targeted DDA with MS1 quantitative methods and analyzing the resulting mass spectrometer data (Adams et al. 2020).
MSstats is an R-based/Bioconductor package for statistical relative quantification of peptides and proteins in mass spectrometry-based proteomic experiments(Choi et al. 2014b). It is applicable to multiple types of sample preparation, including label-free workflows, workflows that use stable isotope labeled reference proteins and peptides, and work-flows that use fractionation. It is applicable to targeted Selected Reactin Monitoring(SRM), Data-Dependent Acquisiton(DDA or shotgun), and Data-Independent Acquisition(DIA or SWATH-MS). This github page is for sharing source and testing.
6.4 Retention Time Correction
For single file, we could get peaks. However, we should make the peaks align across samples for as features and retention time corrections should be performed. The basic idea behind retention time correction is that use the high quality grouped peaks to make a new retention time. You might choose
obiwarp(for dramatic shifts) or loess regression(fast) method to get the corrected retention time for all of the samples. Remember the original retention times might be changed and you might need cross-correct the data. After the correction, you could group the peaks again for a better cross-sample peaks list. However, if you directly use
obiwarp, you don’t have to group peaks before correction.
This paper show a matlab based shift correction methods(Fu et al. 2017). Retention time correction is a Parametric time warping process and this paper is a good start (Wehrens, Bloemberg, and Eilers 2015). Meanwhile, you could use MS2 for retention time correction(Lili Li et al. 2017).
6.5 Filling missing values
Too many zeros or NA in peaks list are problematic for statistics. Then we usually need to integreate the area exsiting a peak.
xcms 3 could use profile matrix to fill the blank. They also have function to impute the NA data by replace missing values with a proportion of the row minimum or random numbers based on the row minimum. It depends on the user to select imputation methods as well as control the minimum fraction of features appeared in single group.
With many groups of samples, you will get another data matrix with column standing for peaks at certain retention time and row standing for samples after the Raw data pretreatment.
6.6 Spectral deconvolution
Without structure information about certain compound, the peak extraction would suffer influence from other compounds. At the same retention time, co-elute compounds might share similar mass. Hard electron ionization methods such as electron impact ionization (EI), APPI suffer this issue. So it would be hard to distinguish the co-elute peaks’ origin and deconvolution method could be used to separate different groups according to the similar chromatogragh behaviors. Another computational tool eRah could be a better solution for the whole process(Domingo-Almenara et al. 2016). Also the ADAD-GC3.0 could also be helpful for such issue(Ni et al. 2016). Other solutions for GC could be found here(Styczynski et al. 2007; Tian et al. 2016b; Du and Zeisel 2013).
6.7 Dynamic Range
Another issue is the Dynamic Range. For metabolomics, peaks could be below the detection limit or over the detection limit. Such Dynamic range issues might raise the loss of information.
Some of the data were limited by the detect of limitation. Thus we need some methods to impute the data if we don’t want to lose information by deleting the NA or 0.
Two major imputation way could be used. The first way is use model-free method such as half the minimum of the values across the data, 0, 1, mean/median across the data(
enviGCMS package could do this via
getimputation function). The second way is use model-based method such as linear model, random forest, KNN, PCA. Try
simputation package for various imputation methods. As mentioned before, you could also use
xcms package to perform imputation.
Tobit regression is preferred for censored data. Also you might choose maximum likelihood estimation(Estimation of mean and standard deviation by MLE. Creating 10 complete samples. Pool the results from 10 individual analyses).
<- rnorm(1000,1) x <0] <- 0 x[x<- x*10+1 y library(AER) <- tobit(y ~ x, left = 0) tfit summary(tfit)
## ## Call: ## tobit(formula = y ~ x, left = 0) ## ## Observations: ## Total Left-censored Uncensored Right-censored ## 1000 0 1000 0 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.0000 0.4398 2.274 0.023 * ## x 10.0000 0.3162 31.623 <2e-16 *** ## Log(scale) 2.1728 0.0000 Inf <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Scale: 8.783 ## ## Gaussian distribution ## Number of Newton-Raphson Iterations: 1 ## Log-likelihood: -3092 on 3 Df ## Wald-statistic: 1000 on 1 Df, p-value: < 2.22e-16
According to Ronald Hites’s simulation(Hites 2019), measurements below the LOD (even missing measurements) with the LOD/2 or with the \(LOD/\sqrt2\) causes little bias and “Any time you have a % non-detected >20%, for whatever reason, it is unlikely that the data set can give useful results.”
Another study find random forest could be the best imputation method for missing at random (MAR), and missing completely at random (MCAR) data. Quantile regression imputation of left-censored data is the best imputation methods for left-censored missing not at random data (Wei et al. 2018).
6.7.2 Over Detection Limit
CorrectOverloadedPeaks could be used to correct the Peaks Exceeding the Detection Limit issue (Lisec et al. 2016).
6.8 RSD/fold change Filter
Some peaks need to be rule out due to high RSD% and small fold changes compared with blank samples.
6.9 Power Analysis Filter
As shown in [Exprimental design(DoE)], the power analysis in metabolomics is ad-hoc since you don’t know too much before you perform the experiment. However, we could perform power analysis after the experiment done. That is, we just rule out the peaks with a lower power for current experimental design.