What are the different marks in the boxplot based on?

In a boxplot, the bottom and top of the box are always the 25th and 75th percentile (the lower and upper quartiles, or Q1 and Q3, respectively), and the band near the middle of the box is always the 50th percentile (the median or Q2). The upper whisker is located at the smaller of the maximum x value and Q3 + 1.5*IQR (Interquantile Range), whereas the lower whisker is located at the larger of the smallest x value and Q1 - 1.5*IQR.

What's the difference between "fold change" and "paired fold change" analyses?

The purpose of fold change (FC) analysis is to compare absolute value change between two group averages. In paired fold change analysis, users aim to find some features that are changed consistently (i.e. up-regulated or down-regulated, but not both) between two groups. The consistency is measured as a percentage - (# of pairs with consistent change above a given FC threshold) / (total # of paired samples). Users need to specify two thresholds - fold change threshold and significant counts (the percentage).

What's the difference between VIP and coefficient - based feature selection (PLS-DA) ?

There are two variable importance measures in PLS-DA. The first, Variable Importance in Projection (VIP), is a weighted sum of squares of the PLS loadings taking into account the amount of explained Y-variation, in each dimension. Please note, since VIP scores are calculated for each components, when more than components are used to calculate the feature importance, the average of the VIP scores are used. The other importance measure is based on the weighted sum of PLS-regression. The weights are a function of the reduction of the sums of squares across the number of PLS components. Please note, for multiple-group (more than two) analysis, the same number of predictors will be built for each group. Therefore, the coefficient of each feature will be different depending on which group you want to predict. The average of the feature coefficients are used as the overall coefficient-based importance.

How to interpret the permutation result (PLS-DA) ?

It is well known that when there are too many variables and a small sample size, many supervised classification algorithms tend to overfit the data. That is, even there is no actual difference between the groups, the program will still be able to discriminate them by picking some features that are "different" between these two group by pure chance! Of course, the classifier will be useless for new data since the pattern it detected is not real (not significant).

The purpose of a permutation test is to answer the question - "what is the performance if the groups are formed randomly". The program uses the same data set with its class labels reassigned randomly. It then builds a new classifer, its performance is then evaluated. The process will repeat many times to estimate the distribution of the performace measure (not necessarily follows a normal distribution). By comparing the performance using the original label and the performance based on the randomly labeled data, one can see if the former is significantly different from the latter. For PLS-DA, the performance is measured using prediction accuracy or group separation distance using the "B/W ratio" (as suggested by Bijlsma et al.). The further away to the right of the distribution formed by randomly permuted data, the more significant the discrimination. The p-value is calculated as the proportion of the times that class separation based on randomly labeled sample is at least as good as the one based on the original data (one-sided p value).

Why the ranks of important features identified by different methods are different ?

Different methods use different criteria for ranking the features. The choice of method used can greatly affect the set of features that are identified. This can be illustrated in the simple cases such as fold change v.s. t-tests, in which the former is interested in absolute change in concentrations/intensities, while the latter focuses on changes relative to the underlying noises. PLS-DA is based on linear regression; SAM and EBAM usually produce similar results since they are all based on t statistics; both random forests and SVM are quite distinctive from all other methods by using an ensemble of classification trees or projections to hyperplanes, respectively. Therefore, these methods tend to generate different results.

How does SAM work?

Significance Analysis of Microarray (SAM) is a robust method designed for identification of statistically significant genes. SAM use moderated t-tests to computes a statistic d_j for each gene j, which measures the strength of the relationship between gene expression (X) and a response variable (Y). This analysis uses non-parametric statistics by repeated permutations of the data to determine if the expression of any gene is significant related to the response. The procedure accounts for correlations in genes and avoids normal assumptions about the distribution of individual genes. More detailed description about about SAM history, features, and instructions can be found here.

Where can I get more information about EBAM?

Empirical Bayesian Analysis of Microarray (EBAM) is based on empirical Bayes method. Both the prior and posterior distributions are estimated from the data. A gene is considered to be differentially expressed if its calculated posterior is larger than or equal to delta and no other genes with a more extreme test score that are not called differentially expressed. The suggested fudge factor a0 is chosen that leads to the largest number of differentially expressed genes. More information about empirical Bayes can be found (here), and detailed description on the EBAM algorithm be found from the original paper by Efron, B., Tibshirani, R. et al.

How does random forest work?

RF is a powerful non-parametric classification method and can be used for both classification and important variable selection (Breiman, 2001). RF uses an ensemble of classification trees, each of which is grown by random feature selection using bootstrap sampling from the original sample set. Class prediction is based on the majority vote of the ensemble. By default, 500 trees are used to build the RF classifier. During tree construction, about one-third of the instances are left out of the bootstrap sampling process. These “left-out” data are then used as test data to obtain an unbiased estimate of the classification error, known as the ‘out-of-bag’ (OOB) error. For a more detailed description about how this algorithm works, click here.

How does recursive SVM (RSVM) work?

This is the workflow of RSVM algorithm described by Zhang X, et al.. Recursive SVM uses SVM for both classification and for selecting a subset of relevant genes according to their relative contribution in the classification. This process is done recursively so that a series of data subsets and classification models can be obtained in a recursive manner, at different levels of feature selection. The performance of the classification can be evaluated either on an independent test data set or by cross validation on the same data set. R-SVM also includes an option for permutation experiments to assess the significance of the performance. Please note, only linear kernel was used for classification, since the information is usually far from sufficient for reliably estimating nonlinear relations for high-dimensional data with a small sample size. First-order approximation will reduce the risk of overfitting in case of limited data. CV2 refers to the cross validation embedded with the feature selection procedure as discussed above.

Can I use the selected features for classification?

Cautions must be taken for the practice. If the features are selected based on the background biological knowledge, it's fine to use the selected subset of data for classification. However, it is not proper to use the features selected based on the whole dataset using some supervised methods (methods that utilize the class label information such as t-tests, SAM, PLS-DA or other supervised classification methods). This procedure will introduce selection bias and result in a very optimistic performance based on cross-validation due to "information leak". In order to obtain an objective performance evaluation, one should include the feature selection procedure in the cross validation. Alternatively, one can evaluate the classifier using independent dataset not used in the feature selections. These procedures were used by PLS-DA, random forest and R-SVM for feature selection and classification. Currently, MetaboAnalyst does not support feature selection and classification using different algorithms (i.e. SAM for feature selection and PLS-DA for classification). This is going to be our next step. Please refer to paper by Ambroise C and McLachlan GJ for a more detailed discussion.

However, one can use non-specific filtering (i.e. not using the class labels) to select features. We recently introduced the Data filter function under the Data Process category. Users can use this function to remove low quanlity features using a variety of criteria. Additionally, users can also try to remove sample outliers if exists using the Data editor. These are safe procedures that can potentially improve the classification performance.

What is a metabolite set?

A metabolite set can be any externally defined groups of metabolites that has something in common. For example,

• Metabolites that are involved in the same pathways (SMPDB, KEGG, Biocarta);
• Metabolites that are changed significantly under the same pathological conditions;
• Metabolites that are associated with the genomic variations (i.e. single nucleotide polymorphism or SNPs);
• Metabolites that are located in the same cellular compartments;
• Other published metabolite signatures, ontology terms, etc.

Please note that you can define your own metabolite set for functional enrichment analysis using MSEA.

What is MSEA?

MESA or Metabolite Set Enrichment Analysis is a novel way to identify biologically meaningful patterns in metabolite concentration changes for quantitative metabolomic studies. In conventional approaches, metabolites are evaluated individually for their significance under conditions of study. Those compounds that have passed certain sigificance level are then combined to see if any meaningful patterns can be discerned. In contrast, MSEA directly investigates if a group of functionally related metabolites are significantly enriched, eliminating the need to preselect compounds based on some arbitrary cut-off threshold. It has the potential to identify "subtle but consistent" changes among a group of related compounds, which may go undetected with conventional approaches.

Essentially, MSEA is a metabolomic version of the popular GSEA (Gene Set Enrichment Analysis) software, with its own collection of metabolite set libraries as well as an implementation of user-friendly web-interfaces. GSEA is widely used in genomics data analysis and has proven to be a powerful alternative to conventional approaches. For more information, please refer to the original paper by Subramanian A, and a nice review paper by Nam D, Kim SY.

How can I obtain more information about these metabolite sets?

Each metabolite set will have an associated source (Database links, PubMed ID, etc) to help users further track down the original source of the information. Please click the "Browse Metabolite Set Libraries" button at the Home page to view all the informations.

The SNP-associated metabolite sets were downloaded from the published journal websites, Only metabolites with p values less than 1e-3 were retained to reduce false positives. For any given SNP ID, users can visit the NCBI SNP database to get detailed information about each SNP.

What are the data input formats supported by enrichment analysis?

There are three types of inputs MSEA accepts:

• A single column data contains a list of compound names for over-representation analysis (ORA). For example, a list of important compounds identified by feature selection or clustering methods;
• A two-column data separated by tab - the first column contains compound name and the second column is their measured concentrations - for single sample profiling (SSP). This data is essentially a single profiled biological sample. The purpose is to compare the measured values to their reference values to see if some compound concentrations are very high or low. This input is limited to human biofluids (Blood, Urine, CSF) for which reference values exist;
• A concentration table. The data contains quantitative metabolomic data of multiple samples measured under different conditions. The data should be uploaded as a comma separated values (.csv) with the each sample per row and each metabolite concentration per column. The first column is sample names and the second column contains sample class labels (discrete or continuous). Users can now directly analyze named metabolite datasets from the Metabolomics Workbench.

What is over-representation analysis (ORA)?

As the name suggests, ORA is to test if certain groups of metabolites are represented more often than expected by chance within a given metabolite list. The most common approach to test this statistically is by using the hypergeometric test or Fisher’s exact test to calculate the probability of seeing at least a particular number of metabolites containing the biological term of interest in the given compound list. MSEA calculates hypergeometrics test score based on cumulative binominal distribution.

What is single sample profiling (SSP)?

One advantage of metabolomics is that metabolite concentrations in biofluids are tightly regulated through homeostasis under normal physiological condition. For common human biofluids such as blood, urine, or cerebral spinal fluids (CSF), normal concentration ranges are known for many metabolites. It is often desirable to know whether certain compound concentrations from a particular sample are significantly higher or lower compared to their normal ranges. For this purpose, we implemented SSP to compare the measured concentrations of compounds to their recorded normal reference range. Compounds that are above or below the normal range beyond a user-specified threshold will then be further investigated using ORA. Please note, SSP is only applicable to human blood, urine and CSF samples.

What is quantitative enrichment analysis (QEA)?

In ORA, only compound names are used for enrichment analysis. In this case, all compounds in the given list are treated equally despite their obvious differences in importance measures (i.e. some are very significant, others only marginally significant, based on p values). This is not optimal. In contrast, QEA takes into consideration not only the counts (presence/absence), but also the importance measures (Q-statistic) calculated for each compounds for enrichment analysis. Therefore, we called quantitative enrichment analysis. In MSEA, QEA was implemented using the globaltest algorithm.

Can I use MSEA on metabolomics data from other species?

Yes, but you may have to provide your self-defined metabolite set libraries. The metabolite set libraries are species specific. For example, the pathway library is only valid for mammalian species. While the metabolite sets based on diseases are only applicable for human. However, MSEA allows users to upload their own metabolite set libraries for enrichment analysis. Therefore, if you want to perform enrichment analysis on data from species other than human being, you need to provide two files, one contains metabolite concentration data, the other contains the metabolite set library.

How to interpret the results from concentration comparisons in single sample profiling (SSP) ?

The reference concentrations SSP were mainly collected from literature. In most cases, the measured concentration should be comparable to the literature reported values as shown below. However, users should keep in mind that these concentrations were measured based on heterogeneous analytical platforms. Therefore, it is more suitable to compare with concentrations measured by the similar technologies. It is advisable to refer to the original literature (provided by SSP) if very extreme values are encountered.

The image below shows the comparisons of the measured urinary Glycine concentration (indicated by a dotted yellow line) to literature reported concentration values. The magenta squares indicate the mean concentration values, and the blue lines indicate the reported ranges. In some cases, 95% confidence intervals are used (2 standard deviations from the mean) if no ranges were provided in the original report. Note the link to the original reference of each study will be provided in a table under each plot during analysis.

What is Q-statistic (QEA) ?

The formula to calculate statistic Q can be obtained from the original publication by (Goeman JJ, et al). Q can be intuitively interpreted as an aggregate of squared covariance between concentration changes and the phenotypes - compounds with large variance have much more influence on the Q than compound with small variance.

How to interpret the metabolite set plot from QEA?

Metabolite set plot shows the detailed influence of individual compound to the total score of the metabolite set. The color of the bar indicates positive or negative association with the phenotype (continuous) or positive association with a particular phenotype (discrete). The reference line at bottom shows the expected influence if the compound was not associated with the phenotype. High bars indicate influential metabolites. The marks on the bars show the standard deviation (SD) of the bar under null hypothesis.

A typical metabolite set plot is shown below. In this case, the phenotype is discrete with four groups. Among all possible comparisons between the data and null hypotheses, the nine amino acids are most significantly associated with either group 3 and 4 as the color indicated. Please note, if you are interested in the comparison between a particular two groups, you should edit and upload your data with samples from only these two groups

Should I use a reference metabolome for my enrichment analysis?

Currently, no single analytical technique can simultaneously measure all the metabolites involved in the metabolic pathways. Different platforms - NMR, GC-MS, LC-MS, usually have different compound coverage. These differences will likely introduce bias during the enrichment analysis, which is not related to the biological questions under investigation. To correct for this bias, a reference metabolome should be used.

A reference metabolome refers to metabolites that can be measured by the analytical platform. Please note, this correction is only necessary when only a list of significant compounds are provided (over-representation analysis). When the whole concentration table is provided, all the compounds within the data are used as the reference metabolome. It is not necessary to upload separately a compound list as reference metabolome.

How is fold-enrichment calculated?

In ORA and SSP: Fold-enrichment is calculated by dividing the observed number of hits by the expected number of hits ("Hits" / "Expect" columns of the ORA and SSP results table).

In QEA: Fold-enrichment is calculated by dividing the observed Q statistic by the expected Q statistic ("Statistic" / "Expected" columns of the QEA results table).

Why aren't all uploaded compound names/IDs matched?

There may be times when not all uploaded compound names/IDs get matches in MetaboAnalyst. This is because our internal compound database only contains all compounds found within our supported pathway/enrichment set libraries.

What is the smart-matching algorithm used to match lipid names?

MetaboAnalyst now implements a "smart-matching" algorithm to help map user's lipid names to our internal compound database. The algorithm considers common lipid synonyms (gathered from PubChem, HMDB and LIPID MAPS), lipid abbreviations (LIPID MAPS), as well as different punctuation/symbols commonly used instead of slashes in lipid names. This step has significantly enhanced the ability of our tool to support lipidomics data for enrichment analysis.

What are the important features for pathay analysis?

There are many commercial pathway analysis software tools, such as Pathway Studio, MetaCore, Ingenuity Pathway Analysis, etc . Compared to them, the pathway analysis in MetaboAnalystis a free, web-based tool targeting for metabolomics data analysis. It uses the high-quality KEGG metabolic pathways as the backend knowledgebase. It integrates many well-established (i.e. univariate analysis, over-representation analysis ) methods, as well as novel algorithms / concepts (GlobalTest, GlobalAncova, pathway topology analysis) into pathway analysis. In addition, MetPA implements a Google-Map style interactive visualization system to help users understand their analysis results.

What is over-representation analysis?

Over-representation analysis is to test if a particular group of compounds is represented more than expected by chance within the user uploaded compound list. In the context of pathway analysis, we are testing if compounds involved in a particular pathway is enriched compared by random hits. The most common methods for such analysis is Fishers' exact test and hypergeometric test. Please note, the over-representation analysis only consider the count (i.e. the total number of compounds that match a particular pathway) and does not consider the magnitude of their concentration changes (not quantitative). So compound that are changed more significant will be treated the same as compounds that are less significant.

What is pathway enrichment analysis?

Pathway enrichment analysis usually refers to quantitative enrichment analysis directly based on the compound concentration values as compared to the compound lists used by over representation analysis. It is usually more sensitive than over-representation analysis and has the potential to discover "subtle but consistent" changes among compounds within the same biological pathway.

Many algorithms have been developed in the last decade for this type of enrichment analysis, the most famous being the Gene Set Enrichment Analysis. Many new and improved methods have since been implemented. The program uses GlobalTest and GlobalAncova for pathway enrichment analysis when users upload concentration tables. Some important features about these two methods include that they support binary, multi-group, as well as continuous phenotypes, and p values can be approximated efficiently based on the asymptotic distribution without using permutations, which is critical for developing web applications.

What is pathway topological analysis?

The structure of biological pathways represents our knowledge about the complex relationships between molecules (activation, inhibition, reaction, etc.). However, neither over-representation analysis or pathway enrichment analysis take the pathway structure into consideration when determining which pathways are more likely to be involved in the conditions under study. It is obvious that changes in the key positions of a network will trigger more severe impact on the pathway than changes on marginal or relatively isolated positions. The program uses two well-established node centrality measures to estimate node importance - betweenness centrality and degree centrality. The former focus on node relative to overall pathway structure, while the latter focus on immediate local connectivities.

Please note, for comparison among different pathways, the node importance values calculated from centrality measures are further normalized by the sum of the importance of the pathway. Therefore, the total/maximum importance of each pathway is 1; the importance measure of each metabolite node is actually the percentage w.r.t the total pathway importance, and the pathway impact is the cumulative percentage from the matched metabolite nodes.

What are the difference between degree centrality and betweenness centrality?

In a graph network, the degree of a node is the number of connections it has to other nodes , Degree centrality is defined as the number of links occurred upon a node. Nodes with higher node degree act as hubs in a network. For directed graph, there are two types of degree: in-degree for links come from other nodes, and out-degree for links initiated from the current node. Metabolic network is directed graph. Here we only consider the out-degree for node importance measure. It is assumed that nodes in upstream will have regulatory roles for the downstream nodes, not vice versa. In the illustration on the right, the nodes colored in red have high degree centrality.

In a graph network, the betweenness centrality measures number of shortest paths going through the node. Therefore, it take into consideration of global network structure, not only immediate neighbour of the current node. For example, nodes that occur between two dense clusters will have a high betweenness centrality, even its degree centrality is not high. Since metabolic network is directed, we use relative betweenness centrality for metabolite importance measure. In the illustration on the right, the nodes colored in blue have high betweenness centrality.

Why should I use a reference metabolome?

Currently, no single analytical technique can simultaneously measure all the metabolites involved in the metabolic pathways. Different platforms - NMR, GC-MS, LC-MS, usually have different compound coverage. These difference will likely introduce bias during the enrichment analysis, which is not related to the biological questions under investigation. To correct for this bias, a reference metabolome should be used. A reference metabolome refers to metabolites that can be measured by the analytical platform.

Please note, this correction is only necessary when only a list of significant compounds are provided (over-representation analysis). When the whole concentration table is provided, the compounds within the data are used as the reference metabolome. It is not necessary to upload separetely a compound list as reference metabolome.

How to interpret the result overview graph (metabolome view)?

An example of "metabolome view" is shown below. It contains all the matched pathways (the metabolome) arranged by p values (from pathway enrichment analysis) on Y-axis, and pathway impact values (from pathway topology analysis) on X-axis. The node color is based on its p value and the node radius is determined based on their pathway impact values. The pathway name is shown as mouse-over tooltip. Click the corresponding node when the tool-tip shows will launch the corresponding pathway view

How to interpret the pathway view?

An example of "pathway view" is shown below. It shows the current metabolic pathway after user clicks the corresponding node on the "metabolome view". The pathway is essentially a simplified KEGG pathway map shown only chemical compounds. The default node color is light blue. When a reference metabolome is provided, some compounds will be colored in light grey if they are not within the reference metabolome. The matched nodes will show varied heat map colors based on their p values. The common compound names can be obtained via mouse over tool-tip. Click on any node will reveal more detailed information as well as database links. For matched compounds, it will also include images summarizing the concentration distribution of the corresponding compound (see compound view).

How to interpret the compound view?

Two examples of "compound view" are shown below. The box plot image on the left shows the result with discrete (binary or multi-group) phenotype labels; the scatter plot image on the right shows the result with continuous phenotype labels. The "compound view" is generated from univariate analysis (t-tests, ANOVA or linear regression). They show the concentration distributions of the corresponding compounds with regard to the phenotype labels.

Why am I unable to view the pathway analysis graphics?

Using pathway analysis, you may find that the graphics no longer appear in your web-browser, but that the results table is still visible. In this case, you may see "Failed to process the request!" appear on the page. To overcome this, delete your web-browser history and cookies, then retry the pathway analysis.

Why aren't all uploaded compound names/IDs matched?

There may be times when not all uploaded compound names/IDs get matches in MetaboAnalyst. This is because our internal compound database only contains all compounds found within our supported pathway/enrichment set libraries.

Why is the number of compounds for a KEGG pathway different between the KEGG website and MetaboAnalyst?

For some KEGG pathways, there may be a discrepancy in the number of compounds listed from KEGG directly compared to MetaboAnalyst. This is because a KEGG pathway can contain compounds that have not been verified for that species. Therefore, our pathway libraries only contain compounds that have been verified so will have fewer compounds than KEGG itself.

How to use and interpret 3D PCA visualization?

The interactive PCA visualization summarizes all the data into the the first three principal components (PCs). Each data point in the Scores Plot represents a sample. Samples that are close together are more similar to each other. The colors of these data points are based on the factor labels. Users can change the colors according to any of the two factor labels.

Each data point in the Loadings Plot represents a feature. When Scores plot and Loadings plot are viewed from the identical perspective, the direction of separation on the scores plot can be explained by the corresponding features on the same directions - i.e. features on the two ends of the direction contribute more to the pattern of separation.

How to use and interpret two-way heatmap?

The two-way heatmap displays the data in the form of matrix with cells colored according to their values. It provides an intuitive overview of the overall data points. To facilitate pattern discovery, MetaboAnalyst also performs clustering on the variables. The samples are sorted by their factor labels.

Users can change the clustering algorithms, the color schema, as well as the order of the samples to suit their preferences. When there are too many samples or variables, the high-resolution image can be used in order to get a clear view of the results.

What is ASCA?

ASCA is a multivariate extension of univariate ANOVA approach. This implementation supports ASCA model for two factors with one interaction effect. The algorithm first partitions the overall data variance (X) into individual variances induced by each factor (A and B), as well as by the interactions (AB). The formula is shown below with (E) indicates the residuals.

X = A + B + AB + E

The SCA part applies PCA to A, B, AB to summarize major variations in each partition. Users need to specify the number of components used for each model. The maximum allowed number of components of each factor must be less than the corresponding levels of the factor. For example, if the phenotype factor contains two levels (control and disease). Only the top component will be extracted. In most cases, users should focus on top one or two components.

How to detect the presence of interaction?

It is critical to test the presence of interactions between the two experimental factors in order to correctly interpret the data.

For individual variable, the interaction can be assessed by univariate two-way ANOVA. Users can view a two-way box plot summary of each variable by clicking the corresponding name in the detailed table view of the two-way ANOVA results.

The overall interaction effect can be assessed by ASCA permutation tests (under the "Model Validation" tab on the ASCA page). It performs unrestricted permutation of observations (Manly's approach) and recalculates the total sum of squares (TSS) for each experimental factor and their interactions. If the TSS calculated based on the original data is significantly different from those calculated from the permuted data, then the effect is significant. For more detailed discussions on this subject, please refer to the paper by Vis D. J. at. al

How to deal with the interaction?

When there are no significant interaction effect being detected, it is advisable to analyze the data with regard to each experimental factor separately for simpler and easier interpretation.

When the presence of interactions between the two experimental factors are confirmed, it is critical to interpret the data in the context of both experimental conditions . Here I will give two examples to illustrate this - one for univariate two-way ANOVA, one for the multivariate ASCA, using the example dataset provided by MetaboAnalyst. The study is to investigate the temporal profiles (factor 1: Time) during a wounding process using two different plant lines (factor 2: Phenotype).

The figure on the left shows the peak (1.4087/417) with very different temporal profiles between the two plant lines. As shown in figure, the baseline (time point 0) of MT is much higher than that of WT group. During a wounding process, WT first shows a significant increase, then starts to decrease; while the MT group shows first decrease, followed by slight increase at the third time point, then decrease again. The different response in the wounding time course may suggest the different underlying biological mechanisms involved.

The figure on the right shows the major trends associated with MT and WT during the wounding time course. The WT first increase then decrease, while the MT first decrease then increase, indicating very different biological processes involved. Please note, the interpretations above are based on the objective of the experiment - to compare temporal profiles of different plant lines. It is also possible to compare the difference between two plant lines at each time point. To summarize, when the interaction is present, researchers must explain the result of factor A in the context of factor B (or vice versa), whichever is more convenient or customary.

How to select (and understand) important variables in ASCA?

ASCA identifies the major trends associated with each experimental factor and their interactions. The next step is to identify variables that are closely follow the detected trends as well as those that clearly diverge from the trends. The first ones would represent those variables that most co-coordinately respond in the experimental context; The second ones would be potential outliers. The variable selection strategy is based on the approach described by Nueda, M.J. et al. An example figure is shown on the left. Important variables are put into two categories: 'well-modeled' refers to those that change co-coordinately to the major profile(s) depicted in the trajectory plots of the factor, and 'outliers' refer to those with relevant changes but different from the major profile(s). These variables are selected based on the Leverage (a measure of the importance of a variable's contribution to the fitted ASCA-model - higher Leverage means more importance) and SPE (or squared prediction error, an evaluation of the goodness of fit of the model to a particular variable - higher SPE means less fit ). The Leverage value (between 0 and 1) is used to select 'well-modeled' variables; while SPE is used to select 'outliers'. SPE is controlled by Alpha (between 0 and 1). Smaller Alpha will lead to higher SPE, which will select variables that diverge more significantly from the general trend.

How to interpret results from multivariate empirical Bayes (MEBA) time-series analysis?

The MEBA approach is designed to compare the time-course profiles with regard to different experimental conditions. It is based on the timecourse method described by YC Tai. et al. The result is a list of variables that are ranked by their difference in temporal profiles across different biological conditions. The Hotelling-T2 is used to rank the variables with different temporal profiles between two biological conditions under study; And the MB-statistics is used for more than two biological conditions. Higher statistical value indicates the time-course profiles are more different across the biological conditions under study.

Can I analyze data with more than two factors?

Although many approaches can be extended to support data with more than two factors. However, as the number of factors increase, the analysis procedures becomes more complicated and data interpretation more difficult. How to present these complicated results in an easy-to-understand form becomes a real challenge for web-based programs (and programmers!). Currently, the task is better left to an expert statistician. We would strongly recommend researchers to further stratify their subjects into more homogenous groups to reduce the number of factors before performing the analysis. In most cases, researchers should focus on one or two factors for better tool support and easier interpretation.

What is the main approach of biomarker analysis offered by MetaboAnalyst?

Receiver Operating Characteristic (ROC) curve analysis is the method of choice for biomarker identification and performance evaluation. However, to perform these two tasks typically requires a good understanding of statistics and machine learning concepts, as well as a good knowledge of programming language such as Matlab or R. MetaboAnalyst offers, through a user-friendly web interface, classical ROC curve analysis, and integration with several well-established algorithms (currently support: PLS-DA, Random Forests and SVM) to assist researchers in performing common ROC curve analysis for biomarker discovery and performance evaluation in metabolomics studies.

What are the recommended sample size?

A minimum of 30-40 samples for balanced groups (i.e. at least 15~20 each group) are recommended in order to calculate decent ROC curves and AUC evaluations. Note, for unbalanced data, the algorithm uses Monte Carlo random sampling to produce balanced sub-samples for training data.

How to understand the ROC curve?

The ROC curve is a fundamental tool for diagnostic test evaluation. It provides a complete and easily visualized sensitivity/specificity report visualization. In a ROC curve the true positive rate (Sensitivity) is plotted in function of the false positive rate (100-Specificity) for different cut-off points of a given parameter. Each point on the ROC curve represents a sensitivity/specificity pair corresponding to a particular decision threshold. A test with perfect discrimination (no overlap in the two distributions) has a ROC curve that passes through the upper left corner (100% sensitivity, 100% specificity). Therefore the closer the ROC curve is to the upper left corner, the higher the overall accuracy of the test. The closer the ROC curve to the diagonal line, the poorer the diagnostic power of the test. The area under the ROC curve is a measure of how well a parameter can distinguish between two diagnostic groups (diseased/normal).

How is the optimal threshold determined in ROC analysis of individual biomarkers?

There are two different approaches to determine the optimal threshold based on ROC curves. In one approach (Youden), the optimal cut-off is the threshold that maximizes the distance to the diagonal line. The optimality criterion is max(sensitivity + specificity). The other approach is to find the point closest to the top-left part of the ROC plot with perfect sensitivity or specificity. The optimality criterion is min((1-sensitivity)^2 + (1-specificity)^2)
Please note, a) there might be more than one threshold identified using the above criteria; b) the optimal cutoff is only calculated for the ROC curve using the actual data points, not the smoothed one.

How are important features selected in multivariate ROC exploratory analysis?

The algorithm tries to identify important features through repeated random sub-sampling cross validation (CV). In each CV, two thirds (2/3) of the samples are used to evaluate the importance of each feature based on VIP scores (PLSDA), decreases in accuracy (Random Forest) or weighted coefficients (Linear SVM). The top 2, 3, 5, 10 ...100 (max) important features are used to build classification/regression models which is validated on the 1/3 the samples that were left out.

Why is is that sometimes the AUROC from the univariate ROC curve analysis is better than those using multiple features?

The ROC curves created (using ROC Explorer or Tester) are based on cross-validation performance of multivariate algorithms (SVM, PLS-DA or Random Forests). In contrast, the classical univariate ROC curves are created based on the performance measured by testing all possible cutoffs within ALL data points. Therefore, the AUROC from cross validated ROC curve is more realistic for prediction purpose, while the AUROC calculated from ROC curve created by univariate approach is often too optimistic. In other words, univariate ROC can be considered as an indicator of the discriminating "potential" of the feature, not its actual performance.

What's the difference between the importance measures in multivariate ROC exploratory analysis?

The most important features are selected from the training data at each cross-validation run. The orders and identities of the list could be (slightly) different each time. There are two forms in reporting the most important features. The Frequency of being selected shows the stability of the rank of the importance for a given biomarker, and the Average importance measure provides a quantitative measure of the importance for a given biomarker. The first measure is more robust and not sensitive to outliers. In most cases, these two approaches should produce exactly the same list.

How do I understand the Prob. Overview?

The prediction overview shows the predicted class probabilities (x-axis) of each sample (y-axis). The probability scores are the average from the 50 iterations, ranging from 0 ~ 1. For instance, less than 0.5 will belong to group A, more than 0.5 belong to group B. In theory, a sample could be located on the 0.5 line, which means the sample has never been selected for testing during the iterations.

Users can also use the Pred. Overview to identify potential outliers. For instance, if a sample is always predicted to have a high probability in the wrong group, this may indicate that the sample could be labeled incorrectly. Users can check "Label samples classified to the wrong groups" to identify these potential outliers.

Why does the result change slightly each time I re-do the analysis?

The algorithm uses repeated random sub-sampling cross validation (CV) to evaluate the feature importance as well as to test the performance of different models. Every time when user click the "Submit" button to re-do the analysis, a new random sampling will be generated. Therefore, the results would be only slightly different. This is the intended behavior, users should looking for the "stable core" for more robust results. The procedure is computationally intensive and the following rules are used to control the time and resources based on the sample size:

• < 100 samples: 50 repeats
• 100 - 200 samples: 30 repeats
• 200 - 500 samples: 20 repeats
• > 500 samples: 10 repeats

Two other factors can also affect the degree of variation - sample size and presence of outliers. If the sample size is small, the traing and testing tend to have high variation. The outliers can also affect the results. Due to the nature of repeated sampling in each CV, some samples may be used multiple times and some samples may never be used. If the outlier samples are used imbalanced (i.e. never been used in the first analysis but used multiple times in the subsequent analysis), the result could change more than "slightly".

How does the program deal with over-fitting?

The algorithm uses repeated random sub-sampling cross validation (CV) to test the performance of models created with different number of features. At each CV, 2/3 samples are used for feature selection and model training and the remaining 1/3 of samples are used for testing. The procedures are repeated 50 times in order to produce a more stable estimation. Users can further perform permutation tests to calculate the significance of the model in ROC Tester. Benchmark testing using different simulated random datasets shows that approach does not overfit. The image below shows some results for classification. The performance is evaluated by measuring the area under ROC curves, error rates, and permutation tests.

How does the program calculate 95% confidence intervals (CIs) of AUC/pAUC?

The program uses non-parametric re-sampling based approach for calculating the confidence intervals. In either bootstrap (classical) or MCCV (multivariate) approach, the data are re-sampled many times, with the AUC/pAUC calculated each time. Confidence intervals are then estimated by simply sorting the data and taking the percentiles to our desired confidence interval bounds. For instance, in 1000 resampling, 95% CIs will be the 25th and 975th values of the sorted AUC values. For details, please refer to the paper by T.A. Lasko et al..

How to understand the permutation tests result?

The permutation tests are performed only for the best model (assuming it containing x features) using the following procedures:

1. Re-assign the phenotype labels randomly to each sample;
2. Perform random sub-sampling cross validation. Within each CV, perform feature ranking and select the top x features to build a classifier using the 2/3 training data, which is then tested on the 1/3 hold-out data. Note, the procedure is repeated only 3 times to save computational time. The performance is then recorded
3. Compare the performance using the original phenotype label and the permuted labels

The performance measures of the permutated data will usually form a normal distribution.If the performance score of the original data lies outside the distribution, then the result is significant. An empirical p value is also usually given. For instance, in 1000 permutation tests, if none of the results are better than the original one, the p value will be reported as p < 0.001.

Note, we have noticed that, when the data does not contain good signals (i.e. AUROC is < 0.65), or when the sample size is small with outliers, the permutation results could be unstable (change from significant to not significant) in different runs. This is unavoidable due to the random subsampling nature of the procedure. Users should always check for performance measures and outliers when this happens.

How to perform the Logistic Regression analysis?

, where α is the intercept term, β is the regression coefficient estimated from the sample dataset, X_i; is the set of covariate (concentration) values, and the P = Pr(y=1|x) is the probability of the outcome of interest (e.g. disease case) as the below.
     P = Pr(Y=outcome of interest | X=x) = exp ( α + βx ) / (1 + exp ( α + βx ) )
The cutoff value of the model will be compared with this predicted probability value if P > cutoff-value, then the sample can be classified as the interest outcome (usually assigned to 1). The best cutoff value will be selected with ROC analysis.

The results of logistic regression including plots and tables are generated using MCCV (100-fold cross-validation) as other methods. In addition, it produces the addition result with 10-fold cross-validation in the "LR model (10-fold CV)" tab.
In order to get the best parsimony model, it may need to try with several combination of compounds/variables, because sometimes it does not work in automatically. However, the best model can be created using the statistics (especially, Lasso Frequency) of variable selection (Builder). The following is the general procedure for building a logistic regression model:

1. Select variables which are most significant using the statistics table (feature ranking) in Builder step.
2. Perform logistic regression analysis with selected variables and get the modeling results.
3. Compare the performance using the accuracy/performance plots and tables (auc, specificity, and sensitivity).
4. Repeat these procedures until you find a satisfied model.

Regarding the class/label of a input data, it is recommended to use the number (0 or 1) instead of string labels. Usually, 1 is used for the case and 0 is for the control.

What is a power analysis?

Power is the probability of detecting an effect, given that the effect is really there. For instance, if a study comparing the two groups (control vs. disease) has power of 0.8 and assuming we can conduct the study many times, then 80% of the time, we would get a statistically significant difference between the two groups, while 20% of time we run this experiment, we will not obtain a statistically significant effect, even though there really is an effect in reality. In practice, researchers are most interested in knowing the sample size (number of subjects) required in order to obtain sufficient power. Power analysis often refers to calculating the sample size required to detect an effect of a given size with a given degree of confidence.

What are the important factors that affect the power?

There are three major factors: (a) effect size which is usually defined as the difference of two group means divided by the pooled standard deviation. When all others are equal, a larger the effect size will lead to more power; (b) degree of confidence which is usually the p value cutoff (alpha) for statistical significance. When all others are equal, there will be reduced power if we require a very high degree of confidence; (c) the sample size - more samples will in general increase power. In many cases, the sample size is our interest for a given power (i.e. 0.8)

How to interpret the diagnostic plots?

The power calculation is based on two assumptions: 1) the effect is indeed present in the data, and 2) the test statistic follow a normal or near normal (Students' t) distribution. The diagnostic plot is to help assess whether these two assumptions are reasonably met. In particular, the panels on the left shows the distribution of the test statistics (t-stat) as a histogram (top panel) and as QQ plot (bottom panel); the panels on the right shows the distribution of the raw p values as histogram (top panel) and the bottom panel shows the sorted p values against their ranks. We expect a large proportion of p values will be close to zero (indicating strong effect) - i.e. the histogram should be left-shifted and the bottom graph should shift to top left corner.

How does MetaboAnalyst compute power?

It is not appropriate to directly apply the traditional univariate approach on power analysis to omics datasets which are characterized by high dimensions with 100s to 1000s of features and a relative small number of samples. Special care needs to be taken in the estimation of effect size as well as the selection of the confidence level. The power analysis in MetaboAnalyst is based on the R package SSPA that was originally developed for power calculation for gene expression datasets. For a given pilot data, the method estimates the average power for a given false discovery rate. For more details, please refer to the original publication by M van Iterson et al

What is the difference between Version 1 and Version 2 of the MS Peaks to Paths module?

With Version 2 of the MS Peaks to Paths module, users can now upload retention time information to perform pathway analysis. Version 2 will only be available 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 Compounds" 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.

How are Empirical Compounds calculated?

Empirical Compounds are intermediaries between m/z features and compounds. The steps for how they are formed are as follows:

1. As in version 1, all m/z features are matched to potential compounds considering different adducts. Then, per compound, all matching m/z features are split into Empirical Compounds based on whether they match within an expected retention time window. The retention time window (in seconds) is calculated as the maximum retention time * 0.02. This results in the initial Empirical Compounds list.
2. Next, Empirical Compounds are merged if they have the same m/z, matched form/ion, and retention time. This results in the merged Empirical Compounds list.
3. Finally, if primary ions are enforced, only Empirical Compounds containing at least 1 primary ion are kept. Primary ions considered are 'M+H[1+]', 'M+Na[1+]', 'M-H2O+H[1+]', 'M-H[-]', 'M-2H[2-]', 'M-H2O-H[-]', 'M+H [1+]', 'M+Na [1+]', 'M-H2O+H [1+]', 'M-H [1-]', 'M-2H [2-]', and 'M-H2O-H [1-]'. This results in the final Empirical Compounds list.
4. Next, pathway libraries are converted from "Compound" space to "Empirical Compound" space. This is done by converting all compounds in each pathway to all Empirical Compound matches. Then the mummichog/GSEA algorithms work as before to calculate pathway enrichment.

What are the inputs for the MS Peaks to Paths module?

Users can upload either peak lists (details below) or a peak intensity table.

Version 1 Peak Lists: The MS Peaks to Pathways module accepts either a three column table containing the m/z features, p-values, and statistical scores, a two-column table containing m/z features and either p-values or t-scores, or a one-column table ranked by either p-values or t-scores. 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.

Version 2 Peak Lists: With Version 2 of the MS Peaks to Pathways module, retention time can be included as a new column with the "rt" or "r.t" heading. The maximum number of columns that can be uploaded is now 5: "m.z", "r.t", "p.value", "t.score" and "mode".

How does the peak intensity table work in the MS Peaks to Paths module?

The peak intensity table is a new addition to the MS Peaks to Paths module that allows for a different way for users to upload their datasets. With the peak intensity table, it is processed (like in the Statistical Analysis module) to create the peak list, and then the mummichog algorithm is performed on that. The only difference is in the visualization of the results. There is now an interactive heatmap that shows the expression of each m/z feature per sample. By default, the m/z features are clustered by p-value. As the visualization is interactive, use your mouse or touchpad to select a part of the heatmap by dragging over a section of the heatmap on the left-side panel (Overview). For instance, this could be a portion that shows an interesting pattern that you wish to further investigate. The selected part will appear in the "Focus View", showing a zoomed-in heatmap. To zoom in or out, change the "Resolution" on the top menu (High to zoom-in, Low to zoom-out). The colors of the heatmap can also be customized using the "Colors" drop-down menu on the top menu. Finally, using the right-hand panel (Enrichment Analysis), perform the mummichog algorithm either on the whole dataset (select entire heatmap), or on a subset of m/z features. Note that if MetaboAnalyst contains multiple libraries for the select organism, the options will appear in the "Database" drop-down menu for users to explore different pathway libraries.

How does the mummichog algorithm for pathway analysis in the MS Peaks to Paths module work?

This module is based on the mummichog algorithm (Li et al.) to perform pathway analysis of untargeted metabolomics data. Briefly, users provide a list of m/z features and a p-value for each feature, highlighting significantly different features between two groups. Three lists are then drawn from this initial list, (1) Lsig, which is the list containing only all significant m/z features (determined by the user selected p-value cutoff), (2) Lref, which is the list of all m/z features, and (3) Lperm, which is a list of randomly drawn m/z features from Lref, but the same length as Lsig. The next steps are as follows:

1. A list of randomly drawn m/z features are drawn from Lref to create Lperm. The m/z features are then mapped to potential metabolites, considering different adducts, protons, etc.
2. The list of potential compounds are then mapped to the user’s selected library of pathways, and a p-value is calculated per pathway.
3. Steps 1 and 2 are repeated many times to compute the null distribution of p-values (modeled as a gamma distribution).
4. The Lsig is mapped to potential metabolites for pathway enrichment analysis, and the resulting p-values (Fisher’s or Hypergeometric, and EASE scores) per pathway for the Lsig compounds are then adjusted for the null-distribution.

Which adducts are used by the MS Peaks to Paths module by default?

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+].

Which currency metabolites are removed by the MS Peaks to Paths module by default?

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

What is the difference between the mummichog algorithm and GSEA?

The original mummichog algorithm implements an over-representation analysis (ORA) method to evaluate pathway-level enrichment based on significant peaks. Users need to specify a pre-defined cutoff based on either t-statistics or fold changes. An alternative approach is the Gene Set Enrichment Analysis (GSEA) method, which is widely used to test for enriched functions from ranked gene lists. Unlike ORA, GSEA considers the overall ranks of features without using a significance cutoff. It can therefore detect subtle and consistent changes which could be missed from using ORA methods.

How are mummichog and GSEA p-values combined?

The MS Peaks to Paths module uses the (Fisher's method) for combining the mummichog and GSEA p-values. It takes the raw p-values per pathway to perform p-value combination.

What are the libraries to select from for pathway analysis?

This module contains both genome-scale metabolic models from the original mummichog python package, KEGG pathway libraries for 21 different organisms, and soon SMPDB pathway libraries for 8 organisms. The genome-scale metabolic model of Danio rerio was manually curated using the KEGG zebrafish model, as well as the human BiGG and Edinburgh Models, and is designated with a [MFN] at the end of its name (Li et al. 2010). The human genome-scale metabolic model was also manually curated and originates from a number of sources (KEGG, BiGG, and Edinburgh Model) and has [MTF] at the end of its name. The remaining four genome-scale metabolic models were directly derived from BioCyc, denoted with [BioCyc]. The KEGG pathway libraries are designated with [KEGG], and the SMPDB pathway libraries will be designated with [SMPDB]. Please select the library closest to your organism.

Why are there several libraries for the same organism?

This module contains both genome-scale metabolic models and KEGG pathway libraries, thereby giving users the power to decide which library to perform their pathway analysis. For instance, for humans, users can either use the manually curated genome-scale metabolic model [MTF], the BioCyc metabolic model [BioCyc], or the KEGG pathway library [KEGG].

Which p-value cut-off should I use? (Mummichog algorithm)

Users have the freedom to select any p-value cut-off to differentiate between significantly and non-significantly enriched m/z features. In general, a cut-off of 0.05 is considered fair, and cut-offs less than 0.05 are considered more stringent. Users should note that their selection of p-value cutoff could affect the results of this module. Changing the p-value cutoff may alter which compounds are included in the Lsig and thereby modify the output of the mummichog algorithm. We therefore recommend that users explore different p-value cutoff values during their analysis, though the algorithm was shown to be robust in the original manuscript (Li et al. 2013).

Why are the gamma adjusted p-values NaN? (Mummichog algorithm)

Users may find that their gamma adjusted p-values are NaN. This is because the permutation p-values could not be Gamma-modeled, and therefore could not be adjusted. This arises because the number of permutations with significance is too close to zero.

How are the p-values adjusted? (Mummichog algorithm)

The results of this module consists of the total number of hits, the raw p-value (Fisher’s), the Gamma-adjusted p-value (for permutations), and the Empirical p-value per pathway. The Gamma-adjusted p-values are the Fisher's p-values calculated using Lsig that are adjusted for the null distribution of permutation p-values. Meanwhile, the Empirical p-values consist of the number of times the p-values using the permuted data were better than the p-values using the significant data, divided by the number of permutations.

What MS instrumentation is supported?

The algorithms used by the MS Peaks to Paths module perform based on the mass accuracy of the instrument (i.e. is not platform specific). However, we recommend that users use high-resolution MS (such as Orbitrap or other Fourier transform mass spectrometry family), which will significantly reduce false positives and improve the pathway activity predictions.

What is the output of this module?

The output of this module first of a plot summarizing the pathway analysis results. Below the plot consists of a table identifying the top-pathways that are enriched in your m/z data. The table is ordered from the most significant to the least significant pathways and contains the number of hits and calculated statistics. A second table is available for download as a link that contains a list of matched metabolites from the user’s uploaded list of m/z features. A third output is the metabolic network visualization, which is based on the KEGG global metabolic network (ko01100).

How do I interpret the metabolic network visualization for this module?

The metabolic network visualization is based on the KEGG global metabolic network (KEGGscape) and has been manually curated. The metabolites are represented as nodes (circles), and their size is based on their associated p-value. More significantly enriched metabolites are bigger then less-significantly enriched metabolites. The aim of this visualization is to provide a global metabolic context for the significantly enriched mapped metabolites, as well as provide an in-house option for visual exploration of the results from the untargeted pathway analysis. Users have the option to colour each specific pathway to their own preference, change the background colour, switch the view style, and download the image as a PNG or SVG file.

What version of the mummichog algorithm is included in MetaboAnalyst?

Currently, MetaboAnalyst supports mummichog version 1.0. In the coming time, we will work to upgrade to mummichog version 2.0, which will use retention-time in the analysis.

How were the networks including metabolite-gene interaction network created?

The networks are generated by first mapping the metabolites, genes or both to the underlying gene-metabolite-disease interactions database. A search algorithm is then performed to identify first-order neighbours (e.g. metabolites that directly interact with a given metabolite) for each of these mapped metabolites ("seeds"). The resulting nodes and their interaction partners are returned to build the subnetworks.

The above figure shows an example of a generated metabolites-genes interaction network.

Which database is used for creating the network?

MetaboAnalyst uses a comprehensive gene-metabolite-disease interaction data based on MetPriCNet. The data contains interaction data mainly extracted from published literature. The database currently contains 11,502 genes, 1,900 metabolites, 132 diseases making a total of 139570, 519, and 78200 interactions for gene-metabolite, metabolite-disease and metabolite-metabolite networks respectively in human. The gene-chemical and chemical-chemical associations for the gene-metabolite and metabolite-metabolite networks, respectively, were extracted from STITCH, such that only highly confident interactions are used. Most of associations in STITCH are based on co-mentions highlighted in PubMed abstracts including reactions from similar chemical structures and similar molecular activities. The associations for the metabolite-disease network were obtained from HMDB. Disease association have been added to HMDB via the Human Metabolome Project’s literature curation team. The Metabolite-Gene-Disease network is an integration of gene-metabolite, metabolite-disease and gene-disease interaction networks.

What is DSPC network and how is it generated?

DSPC stands for debiased sparse partial correlation. It is an algorithm designed to build partial correlation networks based on a graphical LASSO model using a data-driven approach without mapping to any known metabolites interaction databases (Jankova, 2015). DSPC can distinguish between direct and indirect associations and provide insights into the structure of dependencies between metabolites. The main assumption of this DSPC method is that the amount of real metabolites interactions is far smaller than the sample size (i.e., the real partial correlation network among the metabolites is sparse). Under this assumption, DSPC reconstructs a network and calculates partial correlation coefficients and p-values for each pair of metabolic features. Thus, DSPC makes it possible to uncover connectivity patterns among a large number of metabolic features by using fewer samples (Basu et al., 2017). The inferred results can be visualized as weighted networks, where nodes represent the metabolic features, and the edges depict the correlations among them.

The above figure shows an example of a DSPC network, where the red edges represent positive correlations, while the blue edges represent negative correlations

How many nodes can be visualized (network size limit)?

The visualization is actually limited by the performance of users' computers and screen resolutions. Too many nodes will make the network too dense to visualize and the computer slow to respond. We recommend limiting the total number of nodes to between 200 ~ 2000 for the best experience. For very large networks, please make sure you have a decent computer equipped with a modern browser (we recommend the latest Google Chrome).

How do I identify important nodes using degree and betweenness?

Important nodes can be identified based on their position within the network. The assumption is that changes in the key positions of a network will have more impact on the network than changes on marginal or relatively isolated positions. MetaboAnalyst provides two well-established node centrality measures to estimate node importance - degree centrality and betweenness centrality. In a graph network, the degree of a node is the number of connections it has to other nodes. Nodes with higher node degree act as hubs in a network. The betweenness centrality measures the number of shortest paths going through the node. It takes into consideration the global network structure. For example, nodes that occur between two dense clusters will have a high betweenness centrality even if their degree centrality values are not high. Note, you can sort the node table based on either degree or betweenness values by double clicking the corresponding column header.

Can I perform GO or pathway analysis on my selected genes/metabolites alone?

Yes. For query metabolites, you can test for enriched KEGG pathways. For query genes only, you can test enriched gene ontologies or pathways (KEGG/Reactome). To do so, first select and highlight query genes using the Highlight Color toolbar on the top left (you may have to highlight twice for upregulated and downregulated genes respectively); or you can use the Hub Explorer and select queries from the node table. After that, select a functional catergory from the Function Explorer section, and click the Submit button.

How is enrichment analysis performed?

The enrichment analysis is to test whether any functional metabolites/genes from the user selected library are significantly enriched among the currently highlighted nodes within the network. MetaboAnalyst uses hypergeometric tests to compute the enrichment p values.

How do I interpret the differences in node colors and sizes in the default network?

In the default network generated by MetaboAnalyst, the size of the nodes are based on their degree values, with a big size for large degree values. The color of nodes are proportional to their betweenness centrality values. When user switches to Expression View, the color will be based on their expression values (if available).

Can I delete nodes from the network?

Yes. You can delete nodes (with their associated edges) from the current network. First you need to select the nodes from the Node Table in the left pane. Then click the Delete button at the top of the node table. A confirmation dialog will appear asking if you really want to delete these nodes. Note, this action will trigger network re-arrangement, especially if hub nodes are removed. In addition, "orphan" nodes may be produced due to removal. These nodes will also be excluded during re-arrangement.

Can I view all the metabolite/gene members of a pathway within the current graph?

Yes, after you have performed functional enrichment analysis, the over-represented themes will be displayed in the table below. By double clicking on a pathway name, all metabolite/gene members of the pathway will be displayed as highlighted nodes within the current network.

How are the subnetworks generated from my data?

The networks are generated by first mapping the genes/metabolites to the underlying gene-metabolite-disease interaction database. A search algorithm is then performed to identify first-order neighbours (proteins that directly interact with a given protein) for each of these mapped proteins ("seeds"). The resulting nodes and their interaction partners are returned to build the subnetworks. The above approach will typically return one giant subnetwork ("continent") with multiple smaller ones ("islands"). Most subsequent analyses are performed on the continent. Note, networks with less than 3 nodes will be excluded.

What is meta-analysis?

Meta-analysis is a type of statistical technique used to integrate multiple independent datasets that have been collected to study same or similar experimental conditions, in order to obtain more robust biomarkers. By combining multiple data sets, the approach can increase statistical power (more samples) and reduce potential bias.

A key concept in meta-analysis is that it is generally unadvisable to directly combine different independent datasets (i.e. merge them into a single large table) and analyze them as a single unit. This is due to potential batch effects associated with each datasets, which can completely overwhelm the biological effects. This issue has been well-studies in microarray experiments generated from different platforms. It is expected the issue could be more severe due to the lack of standardization in metabolomics.

Instead, meta-analysis is usually computed based on summary statistics (p values, effect sizes, etc.) to identify robust biomarkers. The meta-analysis module in MetaboAnalyst was developed to support these approaches. The results can be visualized using heatmap to explore the patters across different studies.

How to prepare data sets for meta-analysis?

Here are some basic rules for data collection for meta-analysis in MetaboAnalyst:

1. These data sets were collected under comparable experimental conditions, and/or the underlying experiments share the same hypothesis or are held to have the same mechanistic underpinnings;
2. Only two-group comparisons are supported at the moment (i.e., control vs disease);
3. These datasets must share the same type of IDs so that the majority of these feature will match;
4. It is best to keep all data on the same scale or range (i.e. both raw or normalized in the same way). It is generally preferred to compare at log scale. MetaboAnalyst offers boxplots and log normalization to help facilitate this process.

What is p-value combination?

Calculating and combining P-values from multiple studies has long been used in the meta-analysis of microarray data, and we are now using these approaches for metabolomics data. These are simple to calculate and flexible to use. Fisher’s method and Stouffer’s method are two popular approaches. They have similar levels of performance and can be easily interpreted whereby larger scores reflect greater differential abundance. The main difference is that weights (i.e. based on sample sizes) are incorporated into the calculation in Stouffer’s method, whereas Fisher’s method is known as a weight-free method. Note that larger sample size does not warrant larger weights, as the quality of each study can be variable. Users should choose to apply Stouffer’s method only when all studies are of similar quality.

What is vote counting?

Vote counting is the most primitive but simplest and most intuitive method of meta-analysis. In this approach, significant features are first selected based on some criteria (e.g. adjusted P < 0.05 and same direction of fold changes) for each data set. The vote for each features can then be calculated by counting the total number of times it occurs as significant across all data sets. This method is statistically inefficient and should be considered as a last resort in situations when other meta-analysis methods cannot be applied.

When should I use direct merge for meta-analysis?

In this approach, different data sets are merged into a mega-data set and then analyzed as if all data sets were derived from a single experiment. This approach ignores the inherent bias and heterogeneity of data sets from different sources. Many other factors (experiment protocols, technical platforms, raw data processing procedures and so forth) can potentially contribute to the observed differences. This approach should only be used when data sets are very similar (i.e. from the same lab and same platform without batch effects).

Why perform meta-analysis of MS peaks?

Metabolomics is becoming an increasingly popular tool for biological research, however while individual studies may identify certain results, such results may not be reproducible in other independent studies of the same biological questions due to low sample size, sample heterogeneity, the type of LC-MS platform used, or metrics used for interpreting results. Meta-analysis, which is the combination of findings from independent studies - can be used to overcome such limitations and ultimately increase the power, precision, and generalizability of a study.

With the increased use of metabolomics and global efforts towards science transparency, the amount of publicly available metabolomics data deposited in dedicated metabolomics repositories such as Metabolomics Workbench, MetaboLights and OmicsDI has grown tremendously. The potential for researchers in the metabolomics community to enhance their findings with publicly available data is immense, but little effort has been applied for the meta-analysis of untargeted metabolomics data. Such a method would also reduce the bias individual studies may carry towards specific sample processing protocols or LC-MS instruments.

How do you perform meta-analysis of MS peaks?

High-resolution LC-MS has become the dominant platform for global metabolomics. The typical LC-MS workflow starts with data acquisition, followed by data pre-processing, feature selection, compound identification and then pathway analysis. However, peak identification requires significant efforts and is a major challenge for downstream functional interpretation. One strategy is to completely bypass bottleneck of metabolite identification and directly infer functional activity from MS peaks by leveraging the collective knowledge of metabolic pathways and networks. This is the idea of the mummichog algorithm. Essentially what happens is that a list of significant peaks, for instance 169.0814 m/z is matched to the monoistopic mass of potential compounds considering different adducts/isotopes (e.g. Cysteic acid (169.0045) or Norepinephrine (169.0739)).

These putatively matched compounds are then mapped onto known biological pathways. If a list of truly significant features reflect biological activity, then the representation of these true metabolites would be enriched on functional structures (Figure 1), while false matches would be distributed at random (Figure 2). We have extended this idea to support the meta-analysis of untargeted metabolomics data.

Overall, the workflow supports meta-analysis of MS peaks at either a) the pathway level or b) the compound/empirical compound level. Each level of integration has their own merits, for instance at the pathway level, compounds do not have to be matched across all studies, whereas the compound-level integration is better at bringing out weaker signals. The workflow also considers data heterogeneity, meaning that the accuracy of the LC-MS instrument and the MS ionization mode are all taken into account when performing putative metabolite annotation. Ultimately, the goal of this workflow is to enable the reuse/repurposong of multiple datasets to identify a robust meta-signature for the phenotype in question.

One of the bottlenecks for meta-analysis of untargeted metabolomics studies is the heterogeneity in the sample processing, instrumentation used, and even samples themselves that make it difficult to directly compare studies. Another significant hurdle is that unlike targeted metabolomics where a single metabolite can be quantified, in untargeted metabolomics data each metabolite can be linked to multiple peaks, each with their own expression levels. The approach we have adopted is to therefore combine p-values, effect sizes, or both from all peaks across all studies linked to that metabolite to obtain an overall significance level for that metabolite. The workflow supports the use of the original mummichog algorithm, the GSEA algorithm, or both ("integ") to perform the meta-analysis. The meta-analysis can be performed by combining the statistical significance of proposed metabolites before pathway analysis (compound/empirical compound level), or after individual pathway analysis.

What is compound-level meta-analysis of MS peaks?

Using the compound or empirical compound meta-analysis, it will first perform the putative peak annotation on each individual study. Next, the peak annotations are united by default, meaning that only the compounds/empirical compounds found across all studies are kept. Following this step, the list of significant features are updated (see Figure 3 below).

For the mummichog algorithm and "integ" options, a list of significant features are required. Therefore the next step is to combine m/z level statistics across all studies for each compound/empirical compound to obtain the unified list of significant features. Users can either integrate p-values (as if from a T-test), effect sizes, or both to select which m/z features should be used as input (Figure 4). This list is then used as input for pathway activity prediction.

What is pathway-level meta-analysis of MS peaks?

For the pathway-level meta-analysis, each study undergoes the typical steps of calculating m/z level statistics (if not already done), putative metabolite annotation, followed by pathway activity prediction. Once this is complete for all studies, all pathway results are combined to create a unified matrix of pathway-level results (keeping only pathways found across all-studies). Pathway activity scores are then combined using one of several p-value integration methods.

When to use compound or pathway-level meta-analysis?

Pathway-level integration should be used when the studies are independent of eachother - meaning they were performed using different LC-MS instruments, using different extraction methods etc. This is because while lists of features from multiple independent studies looking at the same disease often have little overlap, the use of pathway analysis improves the biological consistency across studies (PMID: 20410053). Towards this, the Pathway-level integration does not require that compounds be matched across all studies. In comparison, the Compound or Empirical-Compound level integration should be used when the data are more homogeneous, for instance data comes from the same lab but different batches or different columns on the same samples and the same instrument. By default, the algorithm will keep only compounds that match across all studies - for instance if the datasets were all the same column, same phenotype in question, same instrument, we would only want to focus on compounds that are merried by all datasets. However, if we had different columns but the same samples, then we would want to keep all putative compounds/empirical compounds for predicting pathway activity.

Which p-value combination method to use?

The meta-analysis of MS peaks framework supports various methods for p-value combination, including Fisher's, Edgington's, Stouffer's, Vote Count, Minimum P-value or Maximum P-value. The choice of statistical method depends on the goal of the meta-analysis as well as the data structure itself. Briefly, Fisher's statistic, the Minimum P-value or Maximum P-value are traditionally used methods for combining p-values.

• The Fisher's method is known to be more sensitive than other methods to very small (or very large) p-values which can result in a high false positive rate. For instance, a single very small p-value can lead to a significant combined p-value. "Fisher’s method employs the log product of individual P-values and thus, a single P-value of zero in one individual case will result in a combined P-value of zero regardless of the other P-values." - (PMID: 26471455). This method should be followed when the data follows a Chi-squared distribution (positive values). It is also not recommended for use for meta-analysis with( >5 datasets).
• The Minimum P-value should be used to answer the question which metabolites are changed across at least one study? In this case, the minimum p-value among all studies is taken as the combined p-value.
• The Maximum P-value is the most restrictive method and should be used to answer the question of which genes are consistently changed across all studies? In this case, the maximum p-value among all studies is taken as the combined p-value.
• Stouffer's method attributes different weights to the p-values when combining them in a meta-analysis. This method should be applied when the data follows a Gaussian curve. This method is not as sensitive as Fisher's to very small or very large p-values.
• Edgington's method uses the sum of the p-values and unlike Fisher's method, is not sensitive to small p-values. It best fits circular data and it has been noted that using this method, a single large p-value can overwhelm small p-values (PMID: 11788962).
• The vote counting method is limited to answering the question is there any evidence of an effect? One issue with this method is that it compares the number of "yes" studies to the number of "no" studies. The "yes" and "no" if often an arbitrary statistical cutoff that can bias the outcome. Secondly, it does not apply any weights to the studies, therefore the effect of a study with 1000 samples has the same weight as a study with 10 samples. While this method is simple to implement and interpret, this method should only be used when standard meta-analysis methods cannot be used. Meanwhile, at least 5 datasets are required for vote counting to reach the statistical confidence. To read more on vote counting: http://assets.press.princeton.edu/chapters/s10045.pdf

Which effect-size combination method to use?

The effect size is a way to quantify the strength of the difference of phenotype in question between two groups. One of the most recommended methods for calculating the effect size is the Hedge's g. This method controls for bias in small studies to overestimate the effect. Cohen's d. To combine the effect sizes across the studies, there are two possible methods: 1) the fixed effects model or 2) the random effects model.

• The fixed effects model is a linear model that assumes that the different studies share an underlying common true effect and that observed differences are due to sampling error. This model should therefore be used when the datasets are rather homogeneous (e.g. samples come from the same population).
• The random effects model assumes that the true effect size follows a distribution, and hence is different between different studies. The random effects model is more commonly used as data are more likely to be heterogeneous (e.g. different LC-MS platform used or batch effects). A fixed effects model will select genes with the strongest effects amongst the studies, while random effects models will select genes with the strongest average effect across the studies.