Inflammatory bowel diseases, which include Crohn’s disease and ulcerative colitis, have affected several million individuals around the world. Jason L. et al. have performed a longtitude multiomics study on the role of microbiome in the pathogenesis of IBD. Metabolomics study on the facal samples is introduced here for example case study of this novel parameters optimization pipeline.
Raw data processing is the first for all metabolomics studies. MetaboAnalystR 3.2 support “.mzXML”, “.mzML” and “.CDF” formats. The original formats (e.g. “.raw” and “.RAW”) generated by vendors will be supported soon. You can convert your original data with ProteoWizard (PMID: 23051804) as the supported formats. Centroided data is preferred.
If you have finished the installation and been ready to use the package. Use the library() function to load the package into R.
# Load the MetaboAnalystR package
library("MetaboAnalystR")
The example Quality Control (QC) samples will be used for the next steps. Download the MS data for parameters' optimization at this step. To reach a goal of quicking learning and avoid the long time running for the whole samples (over 600 samples). We only provide 5 samples for each CD and nonIBD group here. If you want to repeat and verify the results in our manuscript, please go IBD MultiOmics Database, download the full batch data and run them using this pipeline.
## Setting the data depositing folder
data_folder_Sample <- "~/Data_IBD"
data_folder_QC <- "~/QC_IBD"
# Use Google API for data downloading here.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".zip")
# Please authorize your google account to access the data
dl <- drive_download(
as_id("10DBpPEWy2cZyXvmlLOIfpwqQYwFplKYK"), path = temp, overwrite = TRUE)
# Setting your own date file folder
out <- unzip(temp, exdir = data_folder_QC)
# Date files for parameters optimization are deposited below
out
# Now, download the small example data for comparison between CD vs. nonIBD
temp <- tempfile(fileext = ".zip")
dl <- drive_download(
as_id("1-wlFUkzEwWX1afRWLJlY_KEJs7BfsZim"), path = temp, overwrite = TRUE)
# Setting the date file folder
out <- unzip(temp, exdir = data_folder_Sample)
# Date files for normal processing example are deposited below
out
Before running the data analysis, the general data structure and information can be inspected with PerformDataInspect. If there are some extremly significant contaminats, it will be discovered directly.
# Inspect the MS data via a 3D image. "res" are used to specify the resolution for the MS data.
PerformDataInspect(data_folder_QC,res = 50)
## 0020a_XAV_iHMP2_FFA_PREFA02.mzXML
## [1] "RT range is: 18.0985 and 1139.81 seconds !"
## [1] "MZ range is: 69.999908447266 and 849.902038574219 Thomson !"
Instrument specific initial parameters have been embedded. Currently, the most popular platforms that include UPLC-Q/E, UPLC-Q/TOF, UPLC-T/TOF, UPLC-Ion_trap, UPLC-Orbitrap, UPLC-G2S, HPLC-Q/TOF, HPLC-Ion_Trap, HPLC-Orbitrap and HPLC-S/Q. If your platform is not listed here, the “general” option can be applied for this setting. Both “centWave” and “matchedFIlter” modes have been configured for further processing. Besides, the specific parameters could also be set individually, if you have your own opinion about the parameters. This parameters optimization step can be skipped.
# Initial platform specific parameters
param_initial <- SetPeakParam(platform = "UPLC-Q/E")
Parameters optimization with “DoE” strategy. Initial parameters are from the optimized parameters above or the internal default parameters from embedded SetPeakParam function.
# Select relative core according to your work platform
param_optimized <- PerformParamsOptimization(raw_data, param = param_initial, ncore = 8)
# Following results was generated at the Ubuntu 18.04.4
# with 64 GB RAM and Interl Core i7-9700 CPU
# Step 1/6: Start to optimize parameters!
#
# This step may take a long time...
# DoE Optimization Starting Now... (2021-12-20 11:49:19)
# Evaluating Noise level...
# ppm is estimated as : 3.6
# Done!
# Preparing Parameters for optimization finished !
# Data Spliting Finished !
# Peak Preparing Begin...
# Peak Preparing Done !
# Round:1
# DoE Running Begin...
# Round 1 Finished !
# Model Parsing...
# Gaussian peak ratio (%): 51.7
# Model Parsing Done !
#
# Round:2
# DoE Running Begin...
# Round 2 Finished !
# Model Parsing...
# Gaussian peak ratio (%): 54.6
# Model Parsing Done !
#
# Round:3
# DoE Running Begin...
# Round 3 Finished !
# Model Parsing...
# Gaussian peak ratio (%): 50.8
# Model Parsing Done !
#
# No Increase Stopping !
# best parameter settings:
# min_peakwidth: 6.375
# max_peakwidth: 27
# mzdiff: 0.012
# snthresh: 8
# bw: 2
# Peak_method: centWave
# RT_method: loess
# ppm: 3.6
# noise: 6615.83
# prefilter: 2
# value_of_prefilter: 13043.25
# minFraction: 0.8
# minSamples: 1
# maxFeatures: 100
# fitgauss: FALSE
# mzCenterFun: wMean
# integrate: 1
# extra: 1
# span: 0.4
# smooth: loess
# family: gaussian
# verbose.columns: FALSE
# polarity: negative
# perc_fwhm: 0.6
# mz_abs_iso: 0.005
# max_charge: 2
# max_iso: 2
# corr_eic_th: 0.85
# mz_abs_add: 0.001
# rmConts: TRUE
# Step 1/6: Parameters Optimization Finished ! (2021-12-20 12:33:17)
# Time Spent In Total:44mins
The ImportRawMSData function reads in raw MS data files and saves it as an OnDiskMSnExp object. Two plots can be output with SetPlotParam function set with “Plot = T” - the Total Ion Chromatogram (TIC) which provides an overview of the entire spectra, and the Base Peak Chromatogram (BPC) which is a cleaner profile of the spectra based on the most abundant signals. These plots are useful to inform the setting of parameters downstream. For users who wish to view a peak of interest, an Extracted Ion Chromatogram (EIC) can be generated using the PlotEIC function.
# Import raw MS data. The "SetPlotParam" parameters can be used to determine plot or not.
rawData <- ImportRawMSData(NULL, data_folder_Sample, plotSettings = SetPlotParam(Plot=FALSE))
## [1] "The number of CPU cores to be used is set to 6."
## [1] "Successfully imported raw MS data!"
The PerformPeakProfiling function is an updated peak processing pipeline from XCMS R functions that performs peak detection, alignment, and grouping in an automatical step. The function also generates two diagnostic plots including statistics on the total intensity of peaks in different samples, a retention time adjustment map, and a PCA plot showing the overall sample clustering prior to data cleaning and statistical analysis.
# Peak Profiling with optimized parameters
mSet <- PerformPeakProfiling(rawData,param_optimized$best_parameters,
plotSettings = SetPlotParam(Plot = T))
# 6 CPU Threads will be used for peak profiling !
# Step 3/6: Started peak picking! This step will take some time...
# ROI searching under 5 ppm ... Detecting peaks in 73289 regions of interest of 0051_XAV_iHMP2_FFA_SM-9WOBP.mzXML ... OK: 19360 found.
# ROI searching under 5 ppm ... Detecting peaks in 77423 regions of interest of 0393_XAV_iHMP2_FFA_SM-AMR37.mzXML ... OK: 21432 found.
# ROI searching under 5 ppm ... Detecting peaks in 82168 regions of interest of 0407_XAV_iHMP2_FFA_SM-AUP8B.mzXML ... OK: 21737 found.
# ROI searching under 5 ppm ... Detecting peaks in 82894 regions of interest of 0069_XAV_iHMP2_FFA_SM-9OS5Y.mzXML ... OK: 22532 found.
# ROI searching under 5 ppm ... Detecting peaks in 79659 regions of interest of 0433_XAV_iHMP2_FFA_SM-9SNJ4.mzXML ... OK: 24222 found.
# ROI searching under 5 ppm ... Detecting peaks in 75602 regions of interest of 0508_XAV_iHMP2_FFA_SM-9X47O.mzXML ... OK: 19772 found.
# ROI searching under 5 ppm ... Detecting peaks in 77508 regions of interest of 0027_XAV_iHMP2_FFA_SM-77FXR.mzXML ... OK: 24599 found.
# ROI searching under 5 ppm ... Detecting peaks in 75000 regions of interest of 0030_XAV_iHMP2_FFA_SM-6KUCT.mzXML ... OK: 23594 found.
# Step 3/6: Peak picking finished ! (2021-12-20 13:52:57)
#
# Step 4/6: Started peak alignment! This step is running...
# Total of 6239 slices detected for processing... Done !
# Going to the next step...
# PeakGroup -- loess is used for retention time correction.....
# Performing retention time correction using 748 peak groups.
# Applying retention time adjustment to the identified chromatographic peaks ... Done !
# Total of 6239 slices detected for processing... Done !
# Going to the next step...
# Step 4/6: Peak alignment finished ! (2021-12-20 13:55:37)
#
# Step 5/6: Started peak filling! This step may take some time...
# Starting peak filling!
# Defining peak areas for filling-in....OK
# Start integrating peak areas from original files...
# Requesting 636 peaks from 0069_XAV_iHMP2_FFA_SM-9OS5Y.mzXML ... got 337 .
# Requesting 854 peaks from 0393_XAV_iHMP2_FFA_SM-AMR37.mzXML ... got 351 .
# Requesting 849 peaks from 0407_XAV_iHMP2_FFA_SM-AUP8B.mzXML ... got 380 .
# Requesting 1240 peaks from 0051_XAV_iHMP2_FFA_SM-9WOBP.mzXML ... got 412 .
# Requesting 518 peaks from 0027_XAV_iHMP2_FFA_SM-77FXR.mzXML ... got 271 .
# Requesting 716 peaks from 0030_XAV_iHMP2_FFA_SM-6KUCT.mzXML ... got 311 .
# Requesting 711 peaks from 0433_XAV_iHMP2_FFA_SM-9SNJ4.mzXML ... got 309 .
# Requesting 1082 peaks from 0508_XAV_iHMP2_FFA_SM-9X47O.mzXML ... got 365 .
# Step 5/6: Peak filing finished ! (2021-12-20 13:55:58)
# Peak Profiling finished successfully !
# Begin to plotting figures...Done !
Intensity statistics of all samples.
The PerformPeakAnnotation function annotates isotope and adduct peaks with the method from CAMERA (PMID: 22111785). It outputs the result as a CSV file (“annotated_peaklist.csv”) and saves the annotated peaks to mSet object. Finally, the peak list is formatted to the correct structure for MetaboAnalystR and filtered based upon user’s specifications using the FormatPeakList function. This function permits the filtering of adducts (i.e. removal of all adducts except for [M+H]+/[M-H]-) and filtering of isotopes (i.e. removal of all isotopes except for monoisotopic peaks). The goal of filtering peaks is to remove degenerative signals and reduce the file size.
# Setting the Annotation Parameters.
annParams <- SetAnnotationParam(polarity = "negative", mz_abs_add = 0.005)
# Perform peak annotation.
annotPeaks <- PerformPeakAnnotation(mSet, annParams)
# Step 6/6: Starting Peak Annotation...
# Start grouping after retention time.
# Created 520 pseudospectra.
# Generating peak matrix...
# Run isotope peak annotation..
# Found isotopes:584
# Start grouping after correlation...
# Generating EIC's....
# Detecting 0027_XAV_iHMP2_FFA_SM-77FXR.mzXML ... 1187 peaks found!
# Detecting 0030_XAV_iHMP2_FFA_SM-6KUCT.mzXML ... 799 peaks found!
# Detecting 0051_XAV_iHMP2_FFA_SM-9WOBP.mzXML ... 1251 peaks found!
# Detecting 0069_XAV_iHMP2_FFA_SM-9OS5Y.mzXML ... 1421 peaks found!
# Detecting 0393_XAV_iHMP2_FFA_SM-AMR37.mzXML ... 1601 peaks found!
# Detecting 0407_XAV_iHMP2_FFA_SM-AUP8B.mzXML ... 170 peaks found!
# Detecting 0433_XAV_iHMP2_FFA_SM-9SNJ4.mzXML ... 1048 peaks found!
# Detecting 0508_XAV_iHMP2_FFA_SM-9X47O.mzXML ... 674 peaks found!
# Warning: Found NA peaks in selected sample.
# Calculating peak correlations in 520 Groups...
# Calculating graph cross linking in 520 Groups...
# New number of ps-groups: 5785
# mSet has now 5785 groups, instead of 520 !
# Generating peak matrix for peak annotation!
# Polarity is set in annotaParam: negative
# Ruleset could not read from object! Recalculating...
# Calculating possible adducts in 5785 Groups...
# Step 6/6: Successfully performed peak annotation! (2021-12-20 14:03:35)
# Going to the final step...
## Format and filter the peak list for MetaboAnalystR
maPeaks <- FormatPeakList(annotPeaks, annParams, filtIso =F, filtAdducts = FALSE, missPercent = 1)
Part of the table (first 8 rows and 2 samples) are shown below.
| Sample |
X0027_XAV_iHMP2_FFA_SM.77FXR |
X0030_XAV_iHMP2_FFA_SM.6KUCT |
| Label |
CD |
CD |
| 70.4372__894.96 |
891395.64268 |
664966.82422 |
| 70.9412__945.33 |
60450.31055 |
106675.49644 |
| 71.0049__39.16 |
1177400.10095 |
3635910.7839 |
| 72.3855__879.25 |
294771.82782 |
3236.41 |
| 73.0281__52.19 |
2648001.44386 |
266203729.57392 |
| 73.4563__904.89 |
753091.82734 |
165855.85092 |
| 73.9609__947.34 |
133946.14498 |
25463.62699 |
This step can be easily finished with 'Statistical Analysis' module in MetaboAnalyst website. The running R code generated by website could easily be used to repeat the same results. This brief chapter is established to show case how to generate the most results. More statistical utilities can be achieved with the same way.
Remaining in the same working directory, we will read in the filtered peak table. Then, we will replace missing values with a small value (half of the minimum positive value within the original data). Next we will filter varibles out non-informative signals, then normalize the samples to their median and perform log transformation. Finally, we will perform t-test analysis to identify differentially enriched features, setting the FDR-adjusted p-value threshold to 0.25.
# First step is to create the mSet Object, specifying that the data to be uploaded
# is a peak table ("pktable") and that statistical analysis will be performed ("stat").
mSet<-InitDataObjects("pktable", "stat", FALSE)
## Starting Rserve:
## /usr/lib/R/bin/R CMD /home/xialab/R/x86_64-pc-linux-gnu-library/3.6/Rserve/libs//Rserve --no-save
##
##
## R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
## Copyright (C) 2020 The R Foundation for Statistical Computing
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## R is free software and comes with ABSOLUTELY NO WARRANTY.
## You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.
##
## Natural language support but running in an English locale
##
## R is a collaborative project with many contributors.
## Type 'contributors()' for more information and
## 'citation()' on how to cite R or R packages in publications.
##
## Type 'demo()' for some demos, 'help()' for on-line help, or
## 'help.start()' for an HTML browser interface to help.
## Type 'q()' to quit R. ##> SOCK_ERROR: bind error #98(address already in use)
## Rserv started in daemon mode.
## [1] "MetaboAnalyst R objects initialized ..."
# Second step is to read in the filtered peak list, please set the path right first
mSet<-Read.TextData(mSet, "metaboanalyst_input.csv", "colu", "disc")
# The third step is to perform data processing using MetaboAnalystR (filtering/normalization)
# Perform data processing - Data checking
mSet<-SanityCheckData(mSet)
## [1] "Successfully passed sanity check!"
## [2] "Samples are not paired."
## [3] "2 groups were detected in samples."
## [4] "Only English letters, numbers, underscore, hyphen and forward slash (/) are allowed."
## [5] "<font color="orange">Other special characters or punctuations (if any) will be stripped off.</font>"
## [6] "All data values are numeric."
## [7] "A total of 30390 (23.6%) missing values were detected."
## [8] "<u>By default, these values will be replaced by a small value.</u>"
## [9] "Click <b>Skip</b> button if you accept the default practice"
## [10] "Or click <b>Missing value imputation</b> to use other methods"
# Perform data processing - Minimum Value Replacing
mSet<-ReplaceMin(mSet);
# Perform data processing - Variable Filtering and Normalization
mSet<-FilterVariable(mSet, "iqr", "F", 25)
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "MedianNorm", "LogNorm", "NULL", ratio=FALSE, ratioNum=20)
## [1] " Further feature filtering based on Interquantile Range Reduced to 5000 features
## based on Interquantile Range"
## [1] " Further feature filtering based on Interquantile Range Reduced to 5000 features
## based on Interquantile Range"
mSet <- PlotNormSummary(mSet, "norm_0_", "png", 72, width=NA)
mSet <- PlotSampleNormSummary(mSet, "snorm_0_", "png", 72, width=NA)
# The fourth step is to perform fold-change analysis
mSet <- FC.Anal.unpaired(mSet, 2.0, 0)
mSet <- PlotFC(mSet, "fc_0_", "png", 72, width=NA)
# The fifth step is to perform t-test analysis
mSet <- Ttests.Anal(mSet, F, 0.05, FALSE, TRUE)
mSet <- PlotTT(mSet, "tt_0_", "png", 72, width=NA)
## [1] "A total of 218 significant features were found."
# The sixth step is to perform PCA
mSet <- PCA.Anal(mSet)
mSet <- PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet <- PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA, 1,2,0.95,1,0)
mSet <- PlotPCALoading(mSet, "pca_loading_0_", "png", 72, width=NA, 1,2);
mSet <- PlotPCABiplot(mSet, "pca_biplot_0_", "png", 72, width=NA, 1,2)
mSet <- PlotPCA3DScoreImg(mSet, "pca_score3d_0_", "png", 72, width=NA, 1,2,3, 40)
## [1] "The Interactive 3D PCA plot has been created, please find it in mSet$imgSet$pca.3d."
# The seventh step is to perform PLS-DA
mSet <- PLSR.Anal(mSet, reg=TRUE)
mSet <- PlotPLSPairSummary(mSet, "pls_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPLS2DScore(mSet, "pls_score2d_0_", "png", 72, width=NA, 1,2,0.95,1,0)
mSet <- PlotPLS3DScoreImg(mSet, "pls_score3d_0_", "png", 72, width=NA, 1,2,3, 40)
mSet <- PlotPLSLoading(mSet, "pls_loading_0_", "png", 72, width=NA, 1, 2);
mSet <- PLSDA.CV(mSet, "L",5, "Q2")
mSet <- PlotPLS.Classification(mSet, "pls_cv_0_", "png", 72, width=NA)
mSet <- PlotPLS.Imp(mSet, "pls_imp_0_", "png", 72, width=NA, "vip", "Comp. 1", 15,FALSE)
Previous versions of MetaboAnalyst encompassed two modules for functional analysis, metabolic pathway analysis (MetPA) (Xia et al. 2010, https://academic.oup.com/bioinformatics/article/26/18/2342/208464) and metabolite set enrichment analysis (MSEA) (Xia et al. 2010, https://academic.oup.com/nar/article/38/suppl_2/W71/1101310). However, these modules require metabolite identifications prior to use, which remains an important challenge in untargeted metabolomics. In comparison, the mummichog algorithm (Li et al. 2013, http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003123) bypasses the bottleneck of metabolite identification prior to pathway analysis, leveraging a priori pathway and network knowledge to directly infer biological activity based on mass peaks. We have therefore implemented the mummichog algorithm in R in a new module named “MS Peaks to Pathways”. The knowledge-base for this module consists of five genome-scale metabolic models from the original Python implementation which have either been manually curated or downloaded from BioCyc, as well as an expanded library of 21 organisms derived from KEGG metabolic pathways.
To use this module, users must upload a table (.txt) using the Read.PeakListData function containing either:
- Three columns - m/z features, p-values, and t-scores or fold-change values.
- Two columns containing m/z features and either p-values or t-scores
- One column ranked by either p-values or t-scores.
- Four columns - m/z features, p-values, t-scores or fold-change values, and the mode (positive or negative).
All inputted files must be in .txt format. If the input is a three column table, both the mummichog and GSEA algorithms (and their combination) can be applied.
If only p-values (or ranked by p-values) are provided, then only the mummichog algorithm will be applied. If only t-scores (or ranked by t-scores) are provided,
then only the GSEA algorithm will be applied.
If p-values have not yet been calculated, users can use the “Statistical Analysis” module of MetaboAnalyst website or this R package to upload their raw peak tables,
process the data, perform t-tests or fold-change analysis, and then upload these results into the module.
With the table, users also need to specify the type of MS instrument, the ion mode (positive or negative) using the UpdateInstrumentParameters function.
Currently, MetaboAnalystR only supports the handling of peaks obtained from high-resolution MS instruments such as Orbitrap, or Fourier Transform (FT)-MS instruments
as recommended by the original mummichog implementation.
Following data upload, users much select the organism’s library from which to perform untargeted pathway analysis using the SetMass.PathLib function.
Users can then perform the mummichog algorithm on their data using PerformPSEA.
First, users will set algorithm to be used to mummichog using SetPeakEnrichMethod.
Second, users will set the p-value cutoff to delineate between significantly enriched and non-significantly enriched m/z features using the SetMummichogPval function.
Finally, use the PerformPSEA to calcaulte pathway activity.
The output of this module first consists of a table of results identifying the top-pathways that are enriched in the user-uploaded data,
which can be found in your working directory named “mummichog_pathway_enrichment.csv”. The table consists of the total number of hits,
the raw p-value (Fisher’s or Hypergeometric), the EASE score, and the adjusted p-value (for permutations) per pathway.
A second table can also be found in your working directory named “mummichog_matched_compound_all.csv”, that contains all matched metabolites from
the user’s uploaded list of m/z features.
For this tutorial we will first directly use the an example peak list data obtained from untargeted metabolomics of IBD data above (mass accuracy 5.0,
negative ion mode). This is an raw data table generated by the raw data processing pipeline.
Perform T test for further MS peaks to pathway analysis
# Perform t-test
library(MetaboAnalystR);
mSet<-InitDataObjects("pktable", "stat", FALSE);
mSet<-Read.TextData(mSet, "metaboanalyst_input.csv", "colu", "disc")
mSet<-SanityCheckData(mSet)
mSet<-ReplaceMin(mSet);
mSet<-FilterVariable(mSet, "iqr", "F", 25)
feature.nm.vec <- c("")# used to remove certain feature, i.e. c("157.0357")
smpl.nm.vec <- c("")# used to remove certain samples, i.e. c("samples 1")
grp.nm.vec <- c("") # used to remove certain groups, i.e. c("UC")
mSet <- UpdateData(mSet)
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "MedianNorm", "LogNorm", "NULL", ratio=FALSE, ratioNum=20)
mSet<-Ttests.Anal(mSet, F, 0.25, FALSE, TRUE)
# [1] "A total of 1329 significant features were found."
cmSet<-Convert2Mummichog(mSet, rt=TRUE)
mSet<-InitDataObjects("mass_all", "mummichog", FALSE)
SetPeakFormat("mprt")
mSet<-UpdateInstrumentParameters(mSet, 5, "negative");
mSet<-Read.PeakListData(mSet, "mummichog_input_your_date.txt");
mSet<-SanityCheckMummichogData(mSet)
mSet<-SetPeakEnrichMethod(mSet, "integ", version="v2")
# Please set the appropriate p value according to your data
mSet<-SetMummichogPval(mSet, 0.15)
mSet<-PerformPSEA(mSet, "hsa_mfn", "current", permNum =100)
mSet<-PlotPSEAIntegPaths(mSet, "peaks_to_paths_", "png"144)
## [1] "MetaboAnalyst R objects initialized ..."
## [1] "mprt"
## [1] "Retention time tolerance is 22.7888654"
## [1] "Only a small percentage (below 10%) peaks in your input peaks should be significant."
## [2] "The algorithm works best for <u>200~500 significant peaks in 3000~7000 total peaks</u>."
## [3] "A total of 1097 or 21.88 percent signficant mz features were found based on the selected p-value cutoff: 0.2"
## [1] "v1"
## [1] "Got 5014 mass features."
## [1] "1227 matched compounds! cpd2mz"
## [1] "A total of 3636 of duplicates were merged."
## [1] "A total of 3348 of duplicates were merged."
## [1] "A total of 3636 of duplicates were merged."
## [1] "A total of 3636 of duplicates were merged."
## [1] "Resampling, 100 permutations to estimate background ..."
The SetPeakEnrichMethod now has a parameter to let users select whether to use Version 1 or Version 2 of the MS Peaks to Paths algorithms. With Version 2, users can now upload retention time information to perform pathway analysis. Note that Version 2 should only be used if user’s data contains a retention time (“rt” or “r.t”) column. Retention time is used to move pathway analysis from the “Compound” space to “Empirical Compound” space (details below in “How are Empirical Compounds calculated?”). The inclusion of retention time will increase the confidence and robustness of the potential compound matches. Another difference is that currency compounds are removed directly from the user’s selected pathway library, versus removed from potential compound hits during the permutations.
# Perform t-test
mSet<-Ttests.Anal(mSet, F, 0.25, FALSE, TRUE)
## [1] "A total of 1329 significant features were found."
mSet<-Convert2Mummichog(mSet, rt=TRUE)
mSet<-InitDataObjects("mass_all", "mummichog", FALSE)
SetPeakFormat("mprt")
mSet<-UpdateInstrumentParameters(mSet, 5, "negative");
mSet<-Read.PeakListData(mSet, "mummichog_input_your_date.txt");
mSet<-SanityCheckMummichogData(mSet)
mSet<-SetPeakEnrichMethod(mSet, "integ", version="v2")
# Please set the appropriate p value according to your data
mSet<-SetMummichogPval(mSet, 0.15)
mSet<-PerformPSEA(mSet, "hsa_mfn", "current", permNum = 100)
mSet<-PlotPSEAIntegPaths(mSet, "peaks_to_paths_", "png", 144)
##
## [1] "MetaboAnalyst R objects initialized ..."
## [1] "mprt"
## [1] "Retention time tolerance is 22.7888654"
## [1] "Only a small percentage (below 10%) peaks in your input peaks should be significant."
## [2] "The algorithm works best for <u>200~500 significant peaks in 3000~7000 total peaks</u>."
## [3] "A total of 221 or 4.41 percent signficant mz features were found based on the selected p-value cutoff: 0.05"
## [1] "v2" trying URL 'https://www.metaboanalyst.ca/resources/libs/mummichog/hsa_mfn.qs' Content length 634759 bytes (619 KB)
## downloaded 619 KB
##
## [1] "Got 5012 mass features."
## [1] "1226 matched compounds! cpd2mz"
## [1] "2857 inital ECs created!"
## [1] "1430 merged ECs identified!"
##
## Attaching package: igraph
##
## The following objects are masked from package:stats:
##
## decompose, spectrum
##
## The following object is masked from package:base:
##
## union
##
## [1] "352 empirical compounds identified in 1.93101406097412 seconds."
## [1] "A total of 1364 of duplicates were merged."
## [1] "A total of 1269 of duplicates were merged."
## [1] "A total of 1269 of duplicates were merged."
## [1] "A total of 1364 of duplicates were merged."
## [1] "A total of 1087 of duplicates were merged."
## [1] "A total of 1087 of duplicates were merged."
## [1] "A total of 1364 of duplicates were merged."
## [1] "A total of 1269 of duplicates were merged."
## [1] "Resampling, 100 permutations to estimate background ..."
Adduct Customization
Raw MS peaks contain a significant amount of adducts specific to their MS instrument and analytical mode. A comprehensive adduct list is shown below in the “Available” panel. Use this to customize the adduct list used by the Mummichog/GSEA algorithms to match m/z peaks to potential compounds hits. The list of available adducts to choose from can be found in here: positive; here: negative ; here:mixed.
For the negative ion mode, the adducts used are: M-H[-], M-2H[2-], M(C13)-H[-], M(S34)-H[-], M(Cl37)-H[-], M+Na-2H[-], M+K-2H[-], M-H2O-H[-], M+Cl[-], M+Cl37[-], M+Br[-], M+Br81[-], M+ACN-H[-], M+HCOO[-], M+CH3COO[-], and M-H+O[-].
For the positive ion mode, the adducts used are: M[1+], M+H[1+], M+2H[2+], M+3H[3+], M(C13)+H[1+], M(C13)+2H[2+], M(C13)+3H[3+], M(S34)+H[1+], M(Cl37)+H[1+], M+Na[1+], M+H+Na[2+], M+K[1+], M+H2O+H[1+], M-H2O+H[1+], M-H4O2+H[1+], M-NH3+H[1+], M-CO+H[1+], M-CO2+H[1+], M-HCOOH+H[1+], M+HCOONa[1+], M-HCOONa+H[1+], M+NaCl[1+], M-C3H4O2+H[1+], M+HCOOK[1+], and M-HCOOK+H[1+].
Currency Customization
Currency metabolites are abundant substances such as water and carbon dioxide known to occur in normal functioning cells and participate in a large number of metabolic reactions (Huss and Holme, 2007). Because of their ubiquitous nature, removing them can greatly improve pathway analysis and visualization. The list of available currency metabolites can be found here.
By default, the MS Peaks to Paths module considers these metabolites as currency: ‘C00001’, ‘C00080’, ‘C00007’, ‘C00006’, ‘C00005’, ‘C00003’, ‘C00004’, ‘C00002’, ‘C00013’, ‘C00008’, ‘C00009’, ‘C00011’, ‘G11113’, ‘H2O’, ‘H+’, ‘Oxygen’, ‘NADP+’, ‘NADPH’, ‘NAD+’, ‘NADH’, ‘ATP’, ‘Pyrophosphate’, ‘ADP’, ‘Orthophosphate’, and ‘CO2’.
Now we will use another example peak list data obtained from untargeted metabolomics of mice alveolar macrophages in lungs, using an Orbitrap LC-MS (C18 negative ion mode and HILIC positive ion mode). This is an example of the four-column table containing the m.z, mode, p.value, and t.score. We will perform both currency and adduct customization.
# Create objects for storing processed data from the MS peaks to pathways module
mSet<-InitDataObjects("mass_all", "mummichog", FALSE)
# Set peak formart - contains m/z features, p-values and t-scores
SetPeakFormat("mpt")
mSet<-UpdateInstrumentParameters(mSet, 5.0, "mixed");
mSet<-Read.PeakListData(mSet, "https://www.metaboanalyst.ca/MetaboAnalyst/resources/data/mummichog_mixed.txt");
mSet<-SanityCheckMummichogData(mSet)
# Customize currency
curr.vec <-c("Water (C00001)""Proton (C00080)","Oxygen (C00007)","NADPH (C00005)","NADP (C00006)","NADH (C00004)","NAD (C00003)","Adenosine triphosphate (C00002)","Pyrophosphate (C00013)","Phosphate (C00009)","Carbon dioxide (C00011)","Hydrogen (C00282)","Hydrogen peroxide (C00027)","Sodium (C01330)")
# Map selected currency to user's data
mSet<-Setup.MapData(mSet, curr.vec);
mSet<-PerformCurrencyMapping(mSet)
# Now customize adducts
add.vec <-c("M [1+]","M+H [1+]","M+2H [2+]","M+3H [3+]","M+Na [1+]","M+H+Na [2+]","M+K [1+]","M+H2O+H [1+]","M-H2O+H [1+]","M-H4O2+H [1+]","M(C13)+H [1+]","M(C13)+2H [2+]","M(C13)+3H [3+]","M(Cl37)+H [1+]","M-NH3+H [1+]","M-CO+H [1+]","M-CO2+H [1+]","M-HCOOH+H [1+]","M+HCOONa [1+]","M-HCOONa+H [1+]","M+NaCl [1+]","M-C3H4O2+H [1+]","M+HCOOK [1+]","M-HCOOK+H [1+]","M-H [1-]","M-2H [2-]","M-H2O-H [1-]","M-H+O [1-]","M+K-2H [1-]","M+Na-2H [1- ]","M+Cl [1-]","M+Cl37 [1-]","M+HCOO [1-]","M+CH3COO [1-]")
# Set up the selected adducts
mSet<-Setup.AdductData(mSet, add.vec);
mSet<-PerformAdductMapping(mSet, "mixed")
# Perform mummichog algorithm using selected currency and adducts, using Version1 of the mummichog algorithm
mSet<-SetPeakEnrichMethod(mSet, "mum", "v1")
mSet<-SetMummichogPval(mSet, 1.0E-5)
mSet<-PerformPSEA(mSet, "hsa_mfn", "current", 100)
mSet<-PlotPeaks2Paths(mSet, "peaks_to_paths_0_", "png", 72, width=NA) # Image is not shown below to avoid to be overwhelmed.
By default, V2 of the MS Peaks to Pathways algorithms enforces primary ions to be present when creating Empirical Compounds (step 3 above). Users can disable this in the UpdateInstrumentParameters function by setting the “force_primary_ion” parameter to “no”. Additionally, by default the retention-time window used to split Empirical Compounds is calculated as the maximum retention time in the user’s data multiplied by the retention time fraction (default is 0.02). Users can either change the retention time fraction (rt_frac) or set the retention time-window (rt_tol - in seconds). Code details are shown below.
# Disable force primary ion
mSet <-UpdateInstrumentParameters(mSet, instrumentOpt, msModeOpt, force_primary_ion ="no", rt_frac =0.02,rt_tol =NA)
# Change retention time fraction when calculating the retention time window
mSet <-UpdateInstrumentParameters(mSet, instrumentOpt, msModeOpt,force_primary_ion ="yes", rt_frac =0.025, rt_tol =NA)
# Set the retention time window (in seconds)
mSet <-UpdateInstrumentParameters(mSet, instrumentOpt, msModeOpt, force_primary_ion ="yes", rt_frac =0.02, rt_tol =25)
|
Introduction of Batch Effect Correction
The batch effect correction analysis module performs an automatic or specific correction on the peak data table generated by peak picking step. Multiple correction methods have been embedded inside. They include combat, WaveICA, EigenMS, QC_RLSC, ANCOVA, RUV_random, RUV_2, RUVseq series (RUV_s, RUV_r and RUV_g), NOMIS, CCMN. The detailed limitation and mathmatic mechanism of different methods have been illustrated in our manuscript. Here, 4 cases is provided to show users a step to step workflow for both batch effect and signal drift correction.
Data preparation
The peak table generated by FormatPeakList function needs to be manually prepared by hand to supplements the corresponding information below.
-
1st column The samples name;
-
2nd column Sample groups/classes information, if there are any QC samples, please mark them as QC;
-
3th column Batch information, please provide consistent characters within the same batch;
-
4th column order information, please provide injection order information of the whole samples;
-
other column Features’ intensity should be provided. If any internal standards, please mark the feature with “IS”.
NOTE: Please provide as more information as possible. If some information omitted, please leave the columns empty.
Here is an example table.
| S1 |
UC |
B1 |
1 |
3420.067 |
3595.400 |
3403.082 |
3422.175 |
3616.583 |
| S2 |
UC |
B1 |
2 |
3332.439 |
3522.818 |
3666.475 |
3454.279 |
3419.392 |
| S3 |
UC |
B1 |
3 |
3384.851 |
3705.339 |
3620.226 |
3554.620 |
3438.895 |
| S4 |
UC |
B1 |
4 |
3593.815 |
3574.260 |
3620.767 |
3561.398 |
3365.505 |
| QC1 |
QC |
B1 |
5 |
3396.274 |
3407.036 |
3529.572 |
3516.900 |
3383.464 |
| QC2 |
QC |
B1 |
6 |
3400.617 |
3497.153 |
3626.562 |
3573.033 |
3505.492 |
| S5 |
UC |
B1 |
7 |
3507.520 |
3353.179 |
3548.650 |
3428.837 |
3678.415 |
| S6 |
UC |
B1 |
8 |
3454.916 |
3569.033 |
3669.644 |
3312.210 |
3394.256 |
| S7 |
CD |
B1 |
9 |
3628.472 |
3484.903 |
3526.993 |
3655.381 |
3476.667 |
| S8 |
CD |
B2 |
10 |
3639.546 |
3481.560 |
3502.178 |
3724.512 |
3535.429 |
| QC3 |
QC |
B2 |
11 |
3585.290 |
3602.893 |
3546.777 |
3396.528 |
3519.941 |
| QC4 |
QC |
B2 |
12 |
3596.886 |
3474.214 |
3666.306 |
3646.229 |
3473.409 |
| S9 |
UC |
B2 |
13 |
3604.571 |
3594.637 |
3347.312 |
3334.772 |
3419.270 |
| S10 |
CD |
B2 |
14 |
3320.749 |
3421.437 |
3542.846 |
3432.795 |
3431.518 |
This Benchmark data is a large batch data with over 600 samples. The data file size is over 70 MB. This part is designed for repeating our published results. Running of this part may take over 30min to finish. You are strongly recommonded to use your own data instead of this large file to avoid the consuming of patience.
Data Downloading and Preparation
# Use Google API for data downloading peak feature data generated by FormatPeakList here.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".csv")
# Please authorize your google account to access the data
dl1 <- drive_download(
as_id("1wEh2P81J_xFWJs5y4mq98-FsjxJ5wmBO"), path = temp, overwrite = TRUE)
# Use Google API for data downloading meta data here.
# Please "install.packages('googledrive')" and "install.packages('httpuv')"first.
library(googledrive);
temp <- tempfile(fileext = ".csv")
# Please authorize your google account to access the data
dl2 <- drive_download(
as_id("1KaBnSNRrirVPvpRxIubGCqpjX8asNeVA"), path = temp, overwrite = TRUE)
## File downloaded:
## * metaboanalyst_input.csv
## Saved locally as:
## * /tmp/RtmpIBzO2x/file11c64ad429b2.csv
## File downloaded:
## * meta_data.csv
## Saved locally as:
## * /tmp/Rtmp3haVCa/file90a5282e7b0.csv
# Data preparation - read data in & transpose.
# This is a reference example for user to prepare their data.
# Please prepare your data table according to your data format.
MetaboAna_Data <- t(read.csv(dl1$local_path,header = T));
colnames(MetaboAna_Data) <- MetaboAna_Data[1,];
MetaboAna_Data <- MetaboAna_Data[-1,];
MetaboAna_Data <- MetaboAna_Data[order(rownames(MetaboAna_Data)),];
meta_data <- read.csv(dl2$local_path);
meta_data <- meta_data[order(meta_data[,1]),c(1,2,4)];
Prepared_Data <- cbind(meta_data,MetaboAna_Data)[,-4];
write.csv(Prepared_Data,file = "IBD_BC_correction.csv",row.names = F)
datapath <- paste0(getwd(),"/IBD_BC_correction.csv")
Library Package
# Load the MetaboAnalystR package
library("MetaboAnalystR")
Data Filtering and Normalization
mSet<-InitDataObjects("pktable", "stat", FALSE)
## Starting Rserve:
## /usr/lib/R/bin/R CMD /home/qiang/R/x86_64-pc-linux-gnu-library/4.0/Rserve/libs//Rserve --no-save
##
##
## R version 4.0.0 (2020-04-24) -- "Arbor Day"
## Copyright (C) 2020 The R Foundation for Statistical Computing
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## R is free software and comes with ABSOLUTELY NO WARRANTY.
## You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.
##
## Natural language support but running in an English locale
##
## R is a collaborative project with many contributors.
## Type 'contributors()' for more information and
## 'citation()' on how to cite R or R packages in publications.
##
## Type 'demo()' for some demos, 'help()' for on-line help, or
## 'help.start()' for an HTML browser interface to help.
## Type 'q()' to quit R.
##
## Rserv started in daemon mode.
## [1] "MetaboAnalyst R objects initialized ..."
mSet<-Read.TextData(mSet, dl1$local_path, "col", "disc")
mSet<-SanityCheckData(mSet)
mSet<-ReplaceMin(mSet);
mSet<-FilterVariable(mSet, "iqr", "F", 25)
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "MedianNorm", "LogNorm", "NULL", ratio=FALSE, ratioNum=20)
## [1] "Successfully passed sanity check!"
## [2] "Samples are not paired."
## [3] "4 groups were detected in samples."
## [4] "Only English letters, numbers, underscore, hyphen and forward slash (/) are allowed."
## [5] "<font color="orange">Other special characters or punctuations (if any) will be stripped off.</font>"
## [6] "All data values are numeric."
## [7] "A total of 744120 (14.4%) missing values were detected."
## [8] "<u>By default, these values will be replaced by a small value.</u>"
## [9] "Click <b>Skip</b> button if you accept the default practice"
## [10] "Or click <b>Missing value imputation</b> to use other methods"
## [1] " Further feature filtering based on Interquantile Range Reduced to 5000 features based on Interquantile Range"
## [1] " Further feature filtering based on Interquantile Range Reduced to 5000 features based on Interquantile Range"
### Data Orgnization
normalized_set <- mSet[["dataSet"]][["norm"]]
ordered_normalized_set <- normalized_set[order(row.names(normalized_set)), ]
### import metadata
meta_data <- read.csv(dl2$local_path)
new_normalized_set <- cbind(meta_data[,2:4], ordered_normalized_set);
write.csv(new_normalized_set,file = "new_normalized_set.csv")
#perform PCA
mSet <- PCA.Anal(mSet)
mSet <- PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet <- PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet <- PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA,style = 2, 1,2,0.95,0,0)
rm(mSet)
Initializing mSet Object
mSet <- InitDataObjects("pktable", "utils", FALSE)
## [1] "MetaboAnalyst R objects initialized ..."
Data Reading
## we set samples in "row" according to the table format. If your samples are in column, set it as "col".
mSet <- Read.BatchDataTB(mSet, "new_normalized_set.csv", "row")
Automatic Correction
mSet <- PerformBatchCorrection(mSet)
getwd()
## [1] "Correcting using the Automatically method !"
## [1] "Correcting with Combat..."
## [1] "Correcting with WaveICA..."
## [1] "Correcting with EigenMS..."
## [1] "Correcting with QC-RLSC..."
## Cluster size 2216 broken into 551 1665
## Done cluster 551
## Cluster size 1665 broken into 1025 640
## Done cluster 1025
## Done cluster 640
## Done cluster 1665
## [1] "Correcting with ANCOVA..."
## [1] "Best results generated by EigenMS !"
# Show the distances between batches
mSet$dataSet$interbatch_dis
## table combat_edata WaveICA_edata EigenMS_edata QC_RLSC_edata ANCOVA_edata
## 38.49935 39.42686 38.32145 24.18179 36.58138 38.38796
|
