MetabR: an R script for linear model analysis of quantitative metabolomic data
© Ernest et al.; licensee BioMed Central Ltd. 2012
Received: 19 June 2012
Accepted: 8 October 2012
Published: 30 October 2012
Metabolomics is an emerging high-throughput approach to systems biology, but data analysis tools are lacking compared to other systems level disciplines such as transcriptomics and proteomics. Metabolomic data analysis requires a normalization step to remove systematic effects of confounding variables on metabolite measurements. Current tools may not correctly normalize every metabolite when the relationships between each metabolite quantity and fixed-effect confounding variables are different, or for the effects of random-effect confounding variables. Linear mixed models, an established methodology in the microarray literature, offer a standardized and flexible approach for removing the effects of fixed- and random-effect confounding variables from metabolomic data.
Here we present a simple menu-driven program, “MetabR”, designed to aid researchers with no programming background in statistical analysis of metabolomic data. Written in the open-source statistical programming language R, MetabR implements linear mixed models to normalize metabolomic data and analysis of variance (ANOVA) to test treatment differences. MetabR exports normalized data, checks statistical model assumptions, identifies differentially abundant metabolites, and produces output files to help with data interpretation. Example data are provided to illustrate normalization for common confounding variables and to demonstrate the utility of the MetabR program.
We developed MetabR as a simple and user-friendly tool for implementing linear mixed model-based normalization and statistical analysis of targeted metabolomic data, which helps to fill a lack of available data analysis tools in this field. The program, user guide, example data, and any future news or updates related to the program may be found athttp://metabr.r-forge.r-project.org/.
KeywordsR script User-friendly Linear mixed model Statistics Normalization Mass spectrometry-based metabolomics
Quantitative metabolomics is a high-throughput approach to systems biology in which many small molecules (metabolites) from a biological sample are simultaneously measured, commonly using nuclear magnetic resonance spectroscopy (NMR), gas chromatography—mass spectrometry (GC-MS), or liquid chromatography—mass spectrometry (LC-MS). While transcriptomics and proteomics are established approaches for characterizing the effects of experimental conditions on metabolism, gene and protein expression changes merely indicate the potential for changes in metabolic endpoints. Metabolic changes are “real-world” endpoints, so metabolomics can connect these functional genomics platforms with actual physiology.
LC-MS metabolomic approaches fall into two categories: those that attempt to measure every metabolite in the sample (untargeted approaches) and those that attempt to measure only a subset of the metabolites (targeted approaches). A key benefit of targeted approaches is that the detected metabolites can also be readily quantified. Like other approaches to systems biology that rely on the analysis of multiple samples to generate large datasets, two important issues hold true in targeted metabolomics. First, experiments frequently are carried out in multiple “blocks”. For example, targeted LC-MS metabolomic platforms involve lengthy instrumental runs and may rely on multiple runs to enhance metabolite coverage[3, 4], often necessitating multiple run days to analyze all samples. Each run day represents a different block, which introduces technical variability in metabolite detection signals from day-to-day variances in factors related to the instrument’s operation, such as injection volume and ionization efficiency. Second, sampling and measurement variables introduce technical variability in metabolite detection signals, including tissue mass (for multicellular organisms), cell number and size (for microorganisms), sample matrix effects, and mass spectrometer variability (measured by the signal from an internal standard present in the metabolite extraction solvent in our experiments). Clearly, the metabolite signal variability due to block and sampling/measurement variables needs to be distinguished from variability due to experimental treatment factors, which calls for a normalization step to remove the effects of such confounding variables.
Conventional LC-MS metabolomic data normalization is carried out by expressing each metabolite signal relative to values of sampling/measurement variables[3, 4]. Statistical tests for mean differences between treatment groups are performed on normalized metabolite values, with metabolite means averaged across the levels of any block factors (i.e., run day).
There are limitations to this conventional normalization approach, however. First, often many metabolites are normalized to one internal standard (i.e., one for all positive ions and one for all negative ions). This would introduce additional bias if there were low or negative correlation between the internal standard signal and a metabolite signal (i.e., for metabolites with different chemical properties from the internal standard), or if the internal standard signal differed significantly between treatment groups. Second, while ignoring block factors (i.e., comparing metabolite means averaged across samples analyzed on different days) increases sample size, significant block effects on metabolite signals may widen confidence intervals, which may preclude identification of “significant” metabolites and conceal statistical outliers. Block effects may dramatically bias the data, especially if they are not balanced across treatment groups.
Currently available software packages provide powerful tools for pre-processing (i.e., peak selection and integration and retention time alignment), visualization (i.e., biochemical pathway mapping), and/or interpretation of targeted and untargeted metabolomic data[5–10]. However, these packages have limitations because they either 1) do not provide normalization tools for removing confounding effects of experimental variables[7–9]; 2) use the conventional normalization approach; or 3) require the researcher to manually determine a normalization factor for each experimental sample.
A flexible and standardized normalization approach that improves on current limitations would improve metabolomic analyses. An efficient and intuitive approach to control for confounding variables is to estimate their effects on metabolite signals using linear models. Rather than assuming similar relationships between each metabolite signal and confounding variables, a linear model fit for each metabolite can be used to estimate and partition the effects of each experimental variable, including treatment factor, on each metabolite signal. Further, experimental variables can be modeled as having either a fixed or random effect on metabolite signals, with important implications. Fixed-effect variables are assumed to have a constant effect on metabolite signals, influence metabolite signals in an anticipated direction, and have a similar influence in replicate experiments. Common fixed-effect variables are number of cells, tissue mass, and ionization efficiency. By comparison, the effects of random-effect variables cannot be anticipated a priori, and they create variation, but overall do not influence metabolite signals. Typical examples are specimen gender, species or line, experiment day, instrument, and technician, although some of these could be treated as hypothesis-driven experimental factors in some experiments.
Mixed models can be used to estimate the effects of fixed- and random-effect variables on a response variable and are an established approach for normalizing microarray data[12–21]. For two primary reasons, however, currently available microarray data normalization tools are not suitable for metabolomic data. First, microarray normalization tools adjust data for systematic effects specific to microarray technology, such as “dye bias” of different fluors, spatial position effects on the microarray chip, background signals, and biases due to probe binding strengths. Second, microarray normalization tools are often platform specific, designed to carry out pre-processing and quality control only for Illumina BeadArray or for Affymetrix GeneChip platforms, for example.
Given the limitations of current metabolomic data normalization approaches, we developed MetabR, a simple, user-friendly, and stand-alone tool that researchers with no programming background can use to implement linear model-based normalization and statistical analysis of targeted metabolomic data downstream of pre-processing. While MetabR is stand-alone, software with pre-processing tools[5, 6, 8] can be used to generate the input data for MetabR. Further, MetabR generates output files that may be used in subsequent analysis, including normalized data, a heat map and dendrogram, and a comma-separated values (CSV) file formatted for direct upload into Pathway Projector, a web-based biochemical pathway visualization tool.
Implementation of MetabR
where μ = group mean,
Group = treatment factor,
Quantity = a measured, continuous value of the amount of tissue used to produce each sample,
IS = a measured, continuous value of the detection signal from an internal standard present in the metabolite extraction solvent,
Day = a normalization factor accounting for the effects of different run days on metabolite signals,
and e = residual error.
The residuals and treatment group means from the fitted model are added together to yield normalized data, which adjusts for effects of sample quantity, ionization efficiency, and run day, as appropriate for the experimental design of the study.
Output files produced by the MetabR program
Normalized data with technical replicates averaged
A plot of the model residuals for each metabolite vs. each metabolite’s overall mean signal
A plot of the model residuals for each metabolite vs. each metabolite’s overall mean signal, expanded to accommodate metabolite labels
Mean plots for all significant metabolites
Tukey HSD p-values for all treatment group comparisons for every metabolite
q-values for all treatment group comparisons for every metabolite
Mean fold-changes between all treatment group comparisons for every metabolite
Plots of all confounding variables vs. all metabolite measurements, pre- and post-normalization
Heat map and dendrogram of the normalized data
Spreadsheet for direct upload to Pathway Projector
Experimental data collection
Two experimental datasets were generated in our lab to illustrate the utility of MetabR. In both experiments, adipose tissue samples were flash frozen in liquid nitrogen and powdered with a mortar and pestle before metabolite extraction, which followed a previously described procedure. The extracted metabolomes were then analyzed by liquid chromatography—tandem mass spectrometry (LC-MS/MS) via a slightly modified version of the methods of Rabinowitz and co-workers[30–32] that scans for approximately 350 total metabolites in positive and negative ionization modes. The Quan Browser function in the Xcalibur software package (Thermo Scientific, Waltham, MA) was used to assess the presence of each metabolite based on standard detection parameters, such as retention time, signal-to-noise ratio, and peak shape. Signal intensity in ion counts was then determined using Xcalibur to manually integrate each peak, and these data were exported into a Microsoft Excel spreadsheet for statistical analysis.
The first experiment was designed to examine the effects of dietary restriction and insulin immunoneutralization on adipose tissue metabolism in chickens. A total of 127 metabolites were detected in abdominal adipose tissue from 16- or 17-day-old male broiler chicks that were fed ad libitum (“Control”), fasted for 5 hours (“Fast”), or immunoneutralized against the effects of endogenous insulin (“InsNeut”), as we previously described[33, 34]. This study included two factors, Treatment and Day (day 1, day 2, or day 3). Fourteen metabolite measurements from this experiment are provided in Additional files3 (“Chicken example data 1”) and4 (“Chicken example data 2”), corresponding to metabolites detected in positive and negative ionization modes, respectively.
The second experiment was designed to examine the effects of Bisphenol A (BPA) on adipose tissue metabolism in mice. A total of 93 metabolites were detected in abdominal adipose tissue from 32 16-week-old inbred male mice which, from weaning, were fed ad libitum and given drinking water spiked with 0, 0.05, 0.5, or 5 μM BPA. Sixteen mice from each of the inbred strains C57BL/6J and DBA/2J were used in this study. A few missing values arose when a metabolite was not detected in a subset of the samples. Using a zero value for these measurements would bias the results, so they were set to missing (“NA”) which excludes that measurement from analysis. This study included three factors, Treatment, Strain (C57BL/6J or DBA/2J), and Day (day 1, day 2, day 3, or day 4). Twelve metabolite measurements from this experiment are provided in Additional files5 (“Mouse example data 1”) and6 (“Mouse example data 2”), corresponding to metabolites detected in positive and negative ionization modes, respectively.
Modeling confounding variables as fixed- vs. random-effect
In our chicken example, Group, Quantity, and IS were modeled as fixed-effect variables, while Day was modeled as a random-effect variable. To illustrate the difference, if Day is defined as a fixed-effect variable, the estimated treatment group mean includes the average Day effects, and the variance and corresponding confidence intervals are based only on residual error and sample size. Inferences about treatment effects refer only to the days used in the experiment. If Day is defined as a random-effect variable, the estimated mean no longer includes Day. Instead, the Day effect becomes a source of random variation that is added to the variance of the estimated mean. The variance and confidence intervals are larger than those when Day is treated as a fixed-effect variable, but experimental results can now be correctly extrapolated to all possible days.
Chicken experimental results
For the chicken data, Quantity (tissue mass) and IS (internal standard measurement, Tris in positive ionization mode and Benzoic Acid in negative ionization mode) were selected as fixed-effect regression variables, and Day (run day) as a random-effect factor.
Chicken experiment fold-changes
Table 2 contains all between-group mean fold-changes for the metabolites, with differences tested by Tukey’s HSD at the 5% significance level. We produced this table by combining the mean fold-changes and p-values exported automatically by MetabR. As shown, the experiment had sufficient power to detect a fold-change as low as 1.25 for Citrate between Fast and Control groups. In general, the differences between the Control and InsNeut groups were smaller than other treatment group comparisons. The program exports q-values automatically, and the researcher may select p-value, q-value, mean fold-change, or a combination of either p-value or q-value and mean fold-change as a significance threshold. As technological improvements continue to allow more metabolites to be detected, the chance of false discoveries will increase, making false discovery corrections (q-value) increasingly necessary.
Mouse experimental results
Mouse experiment fold-changes
The open-source statistical computing software R provides a convenient environment for statistical analysis of metabolomic and other -omic data. We developed a user-friendly R program that normalizes metabolomic data using linear mixed-effect modeling (with regression and ANOVA terms), statistically compares treatments, and exports results files to aid data interpretation, filling an important lack in statistical analysis tools available to the metabolomics community. The MetabR program file, example data, and user guide are available as an R-Forge project athttp://metabr.r-forge.r-project.org/. This website will also contain future news or updates related to MetabR, including availability through the Comprehensive R Archive Network (CRAN) or Bioconductor.
Availability and requirements
Project name: MetabR
Project home page: http://metabr.r-forge.r-project.org/
Operating system(s): Windows, Mac, Linux, any system that runs R
Programming language: R
Other requirements: Required R packages are installed automatically. The program was written and tested using R version 2.15 for Windows.
License: GNU General Public License (GPL)
Any restrictions to use by non-academics: No restrictions
Availability of supporting data
The datasets supporting the results of this article are included within the article (and its additional files).
J Gooding’s current address: Sarah W. Stedman Nutrition & Metabolism Center, Duke University School of Medicine, 4321 Medical Park Drive, Suite 200, Durham, NC 27704
Analysis of variance
Graphical user interface
Honest Significant Difference
liquid chromatography—mass spectrometry
Liquid chromatography—tandem mass spectrometry.
JRG and SRC were supported by funding from the National Science Foundation through an Ocean Sciences award (OCE-1061352) to the University of Tennessee at Knoxville. Funding for metabolomic analyses of chicken adipose tissue was provided by a University of Tennessee AgResearch Innovation Grant to BHV and SRC.
The authors thank Brantley Wyatt, previously of the University of Tennessee Graduate School of Genome Science and Technology, for conducting the mouse experiments and generating the mouse adipose tissue samples used in this work, and Drs. Joelle Dupont and Jean Simon of the Institut National de la Recherche Agronomique (INRA) for conducting the chicken experiments and providing the corresponding adipose tissue samples.
- Nicholson JK, Connelly J, Lindon JC, Holmes E: Metabonomics: a platform for studying drug toxicity and gene function. Nat Rev Drug Discov. 2002, 1: 153-161. 10.1038/nrd728.PubMedView ArticleGoogle Scholar
- Reaves ML, Rabinowitz JD: Metabolomics in systems microbiology. Curr Opin Biotechnol. 2011, 22: 17-25. 10.1016/j.copbio.2010.10.001.PubMedPubMed CentralView ArticleGoogle Scholar
- Tai E, Tan M, Stevens R, Low Y, Muehlbauer M, Goh D, Ilkayeva O, Wenner B, Bain J, Lee J, Lim S, Khoo C, Shah S, Newgard C: Insulin resistance is associated with a metabolic profile of altered protein metabolism in Chinese and Asian-Indian men. Diabetologia. 2010, 53: 757-767. 10.1007/s00125-009-1637-8.PubMedPubMed CentralView ArticleGoogle Scholar
- Kwon YKI, Higgins MB, Rabinowitz JD: Antifolate-induced depletion of intracellular glycine and purines inhibits thymineless death in E. coli. ACS Chem Biol. 2010, 5: 787-795. 10.1021/cb100096f.PubMedPubMed CentralView ArticleGoogle Scholar
- Xia J, Psychogios N, Young N, Wishart DS: MetaboAnalyst: a web server for metabolomic data analysis and interpretation. Nucleic Acids Res. 2009, 37: W652-W660. 10.1093/nar/gkp356.PubMedPubMed CentralView ArticleGoogle Scholar
- Creek DJ, Jankevics A, Burgess KEV, Breitling R, Barrett MP: IDEOM: an Excel interface for analysis of LC–MS-based metabolomics data. Bioinformatics. 2012, 28: 1048-1049. 10.1093/bioinformatics/bts069.PubMedView ArticleGoogle Scholar
- Smith CA, Want EJ, O’Maille G, Abagyan R, Siuzdak G: XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification. Anal Chem. 2006, 78: 779-787. 10.1021/ac051437y.PubMedView ArticleGoogle Scholar
- Melamud E, Vastag L, Rabinowitz JD: Metabolomic analysis and visualization engine for LC − MS data. Anal Chem. 2010, 82: 9818-9826. 10.1021/ac1021166.PubMedView ArticleGoogle Scholar
- Kono N, Arakawa K, Ogawa R, Kido N, Oshita K, Ikegami K, Tamaki S, Tomita M: Pathway Projector: web-based zoomable pathway browser using KEGG Atlas and Google Maps API. PLoS One. 2009, 4: e7710-10.1371/journal.pone.0007710.PubMedPubMed CentralView ArticleGoogle Scholar
- Boccard J, Veuthey JL, Rudaz S: Knowledge discovery in metabolomics: an overview of MS data handling. J Sep Sci. 2010, 33: 290-304. 10.1002/jssc.200900609.PubMedView ArticleGoogle Scholar
- Oberg L, Mahoney DH: Linear mixed effects models. Topics in Biostatistics. Edited by: Ambrosius WT. 2007, Totowa, NJ: Humana Press, 213-234.View ArticleGoogle Scholar
- Wolfinger RD, Gibson G: Assessing gene significance from cDNA microarray expression data via mixed models. J Comput Biol. 2001, 8: 625-637. 10.1089/106652701753307520.PubMedView ArticleGoogle Scholar
- Yang YH, Dudoit S, Luu P, Speed TP: Normalization for cDNA microarry data. SPIE Proceedings. 2001, 4266: 141-152.View ArticleGoogle Scholar
- Berger MPF, Passos VL, Tan FES, Winkens B: Optimal designs for one- and two-color microarrays using mixed models: a comparative evaluation of their efficiencies. J Comput Biol. 2009, 16: 67-83. 10.1089/cmb.2008.0048.PubMedView ArticleGoogle Scholar
- Chu T-M, Weir B, Weir , Wolfinger R: A systematic statistical linear modeling approach to oligonucleotide array experiments. Math Biosci. 2002, 176: 35-51. 10.1016/S0025-5564(01)00107-9.PubMedView ArticleGoogle Scholar
- Demirkale CY, Nettleton D, Maiti T: Linear mixed model selection for false discovery rate control in microarray data analysis. Biometrics. 2010, 66: 621-629. 10.1111/j.1541-0420.2009.01286.x.PubMedView ArticleGoogle Scholar
- Haldermans P, Shkedy Z, Van Sanden S, Burzykowski T, Aerts M: Using linear mixed models for normalization of cDNA microarrays. Stat Appl Genet Mol Biol. 2007, 6:
- Li H, Wood C, Getchell T, Getchell M, Stromberg A: Analysis of oligonucleotide array experiments with repeated measures using mixed models. BMC Bioinforma. 2004, 5: 209-10.1186/1471-2105-5-209.View ArticleGoogle Scholar
- Wang L, Zhang B, Wolfinger RD, Chen X: An integrated approach for the analysis of biological pathways using mixed models. PLoS Genetics. 2008, 4: e1000115-10.1371/journal.pgen.1000115.PubMedPubMed CentralView ArticleGoogle Scholar
- Urs S, Smith C, Campbell B, Saxton AM, Taylor J, Zhang B, Snoddy J, Jones Voy B, Moustaid-Moussa N: Gene expression profiling in human preadipocytes and adipocytes by microarray analysis. J Nutr. 2004, 134: 762-770.PubMedGoogle Scholar
- Wernisch L, Kendall SL, Soneji S, Wietzorrek A, Parish T, Hinds J, Butcher PD, Stoker NG: Analysis of whole-genome microarray replicates using mixed models. Bioinformatics. 2003, 19: 53-61. 10.1093/bioinformatics/19.1.53.PubMedView ArticleGoogle Scholar
- Smyth GK, Speed T: Normalization of cDNA microarray data. Methods. 2003, 31: 265-273. 10.1016/S1046-2023(03)00155-5.PubMedView ArticleGoogle Scholar
- Du P, Kibbe WA, Lin SM: lumi: a pipeline for processing Illumina microarray. Bioinformatics. 2008, 24: 1547-1548. 10.1093/bioinformatics/btn224.PubMedView ArticleGoogle Scholar
- Verzani J: An introduction to gWidgets. R News. 2007, 7: 26-33.Google Scholar
- Bates D, Maechler M, Bolker B: lme4: Linear mixed-effects models using S4 classes. 2011, [http://CRAN.R-project.org/package=lme4],Google Scholar
- R Development Core Team: R: A language and environment for statistical computing. R Foundation for Statistical Computing. 2011, Vienna, Austria: R Foundation for Statistical Computing, [http://www.R-project.org/],Google Scholar
- Noguchi K, Hui WW, Gel YR, Gastwirth JL, Miao W: lawstat: An R package for biostatistics, public policy, and law. 2009, [http://CRAN.R-project.org/package=lawstat],Google Scholar
- Storey JD: A Direct approach to false discovery rates. Journal of the Royal Statistical Society B. 2002, 64: 479-498. 10.1111/1467-9868.00346.View ArticleGoogle Scholar
- Saxton AM: A macro for converting mean separation output to letter groupings in Proc Mixed. 1996, Nashville: Proceedings, 23rd SAS Users Group International: 22-25 March 1998, 1243-1246.Google Scholar
- Collier JJ, Burke SJ, Eisenhauer ME, Lu D, Sapp RC, Frydman CJ, Campagna SR: Pancreatic β-cell death in response to pro-inflammatory cytokines Is distinct from genuine apoptosis. PLoS One. 2011, 6: e22485-10.1371/journal.pone.0022485.PubMedPubMed CentralView ArticleGoogle Scholar
- Bajad SU, Lu W, Kimball EH, Yuan J, Peterson C, Rabinowitz JD: Separation and quantitation of water soluble cellular metabolites by hydrophilic interaction chromatography-tandem mass spectrometry. Journal of Chromatography A. 2006, 1125: 76-88. 10.1016/j.chroma.2006.05.019.PubMedView ArticleGoogle Scholar
- Waters CM, Lu W, Rabinowitz JD, Bassler BL: Quorum sensing controls biofilm formation in Vibrio cholerae through modulation of cyclic di-GMP levels and repression of vpsT. J Bacteriol. 2008, 190: 2527-2536. 10.1128/JB.01756-07.PubMedPubMed CentralView ArticleGoogle Scholar
- Dupont J, Tesseraud S, Derouet M, Collin A, Rideau N, Crochet S, Godet E, Cailleau-Audouin E, Metayer-Coustard S, Duclos MJ, Gespach C, Porter TE, Cogburn LA, Simon J: Insulin immuno-neutralization in chicken: effects on insulin signaling and gene expression in liver and muscle. J Endocrinol. 2008, 197: 531-542. 10.1677/JOE-08-0055.PubMedView ArticleGoogle Scholar
- Ji B, Ernest B, Gooding J, Das S, Saxton A, Simon J, Dupont J, Metayer-Coustard S, Campagna S, Voy B: Transcriptomic and metabolomic profiling of chicken adipose tissue in response to insulin neutralization and fasting. BMC Genomics. 2012, 13: 441-10.1186/1471-2164-13-441.PubMedPubMed CentralView ArticleGoogle Scholar
- Warnes GR: gplots: Various R programming tools for plotting data. 2012, [http://CRAN.R-project.org/package=gplots],Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.