Functional Analysis of Global Metabolomics

2.1 Peak list

To use this module, users have multiple options to upload their data. The typical data formats is a table (*.txt) containing one to four columns. Users can format the table based on the following rules,

1). Three columns - m/z features, p-values, and t-scores or fold-change values.

2). Two columns containing m/z features and either p-values or t-scores

3). One column ranked by either p-values or t-scores.

4). Four columns - m/z features, p-values, t-scores or fold-change values, and the mode (positive or negative).

Data uploading can be performed by using the Read.PeakListData function.

If p-values have not yet been calculated, users can use the “Statistical Analysis” module from MetaboAnalyst website or MetaboAnalystR (see more details from Statistical Analysis (one-factor) vignette) to upload their raw peak tables, process the data, perform t-tests, 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), and the p-value cutoff to delineate between significantly enriched and non-significantly enriched m/z features using the UpdateMummichogParameters 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 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 calculate pathway activity.

See an example from section 3 (Case Study 1) in this vignette.

2.2 Peak table

MetaboAnalystR also accepts a typical peak intensity table, including complete features’ information, groups of samples and peak height(or areas). Here is the basic rules for peak table:

1). Feature information must contains both m/z and retention time (RT). They should be pasted a one string with double underscore (__) . e.g. m/z 157.0241 and RT 28.64 sec should be formatted as “157.0241__28.64”. See example below.

2). At least two groups are required for functional analysis. User may provide more groups, but it is highly recommended to include only TWO groups each time to discover the functional perturbation;

3). Intensity values must be a number, all NAs will be replaced with 0;

Here is an example peak intensity table.

Different from peak list option, there are a few more steps for peak intensity feature processing. In details, Read.TextData is used to read the peak table. ReplaceMin and FilterVariable function is used to replace the minimum value and filter out the redundancy. PreparePrenormData and Normalization is also needed to normalize your peak table. Finally, execute PreparePeakTable4PSEA function to automatically convert peak table into compatible peak list for the following processing. Other steps are similar to the ones implemented for peak list.

See an example from section 4 (Case Study 2) in this vignette.

2.3 MS/MS identification results

Since the version 4.0, MetaboAnalystR has also accepted the MS/MS based identification results, together with MS peak list for more accurate pathway prediction. To prepare the MS peak and MS/MS based identification data, user need to follow some rules to format their data. Here are two files needed for MS/MS identification based functional analysis.

1). MS peak list is same as the one described in section 2.1. But users are highly encouraged to include four columns (mz, rt, t.score and p.value);

2). MS/MS compound list must include same number of rows as the MS peak list. The header of each column can be any characters based on users’ preference. Besides, the compound identifications can be one of these types (InChiKeys, PubChem CIDs, PubChem SIDs and SMILES). Other compound ID will be supported soon;

To use MS/MS based identification results, some extra functions need to be executed together with the typical mummichog functions. In details, SetMS2IDType is used to set ID type of compounds from MS/MS identification. Different from section 2.1 and 2.2, Read.PeakMS2ListData is used to import both MS features and MS/MS identification results consurrently. Other functions and steps are similar to the ones used to process MS peaks in section 2.1.

See an example from section 5 (Case Study 3) in this vignette.

2.4 Mummichog version 1 or version 2?

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 (see details 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.

2.5 Results output

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_mummichog.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 directly use several examples collected using LC-high resolution MS) to showcase the performance and implementation.

6.1 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 positive; negative; 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+].

6.2 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.

# Load MetaboAnalystR
library(MetaboAnalystR)

# Clean global environment
rm(list = ls())

# Create objects for storing processed data
mSet <- InitDataObjects("mass_all", "mummichog", FALSE)

## [1] "MetaboAnalyst R objects initialized ..."

# Set peak formart - contains m/z features, p-values and t-scores
mSet <- SetPeakFormat(mSet, "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);

## [1] "Loaded files from MetaboAnalyst web-server."

# 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");

## [1] "Loaded files from MetaboAnalyst web-server."

# 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);

## [1] "Got 3915 mass features."
## [1] "A total of 1188 of duplicates were merged."
## [1] "A total of 1292 of duplicates were merged."
## [1] "A total of 1188 of duplicates were merged."
## [1] "A total of 1188 of duplicates were merged."
## [1] "Resampling,  100 permutations to estimate background ..."

# Image is not shown below to avoid to be overwhelmed.
mSet <- PlotPeaks2Paths(mSet, "peaks_to_paths_2_", "png", 72, width=NA)