Abstract
Background
Metabolomics is an emerging highthroughput 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 fixedeffect confounding variables are different, or for the effects of randomeffect 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 randomeffect confounding variables from metabolomic data.
Findings
Here we present a simple menudriven program, “MetabR”, designed to aid researchers with no programming background in statistical analysis of metabolomic data. Written in the opensource 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.
Conclusions
We developed MetabR as a simple and userfriendly tool for implementing linear mixed modelbased 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 at http://metabr.rforge.rproject.org/ webcite.
Keywords:
R script; Userfriendly; Linear mixed model; Statistics; Normalization; Mass spectrometrybased metabolomicsFindings
Background
Quantitative metabolomics is a highthroughput 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 (GCMS), or liquid chromatography—mass spectrometry (LCMS). 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 “realworld” endpoints, so metabolomics can connect these functional genomics platforms with actual physiology [1].
LCMS 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) [2]. 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 LCMS 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 daytoday 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 LCMS 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 preprocessing (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 [510]. However, these packages have limitations because they either 1) do not provide normalization tools for removing confounding effects of experimental variables [79]; 2) use the conventional normalization approach [6]; or 3) require the researcher to manually determine a normalization factor for each experimental sample [5].
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. Fixedeffect 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 fixedeffect variables are number of cells, tissue mass, and ionization efficiency. By comparison, the effects of randomeffect 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 [11], although some of these could be treated as hypothesisdriven experimental factors in some experiments.
Mixed models can be used to estimate the effects of fixed and randomeffect variables on a response variable [11] and are an established approach for normalizing microarray data [1221]. 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 [22]. Second, microarray normalization tools are often platform specific, designed to carry out preprocessing and quality control only for Illumina BeadArray or for Affymetrix GeneChip platforms, for example [23].
Given the limitations of current metabolomic data normalization approaches, we developed MetabR, a simple, userfriendly, and standalone tool that researchers with no programming background can use to implement linear modelbased normalization and statistical analysis of targeted metabolomic data downstream of preprocessing. While MetabR is standalone, software with preprocessing 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 commaseparated values (CSV) file formatted for direct upload into Pathway Projector [9], a webbased biochemical pathway visualization tool.
Methods
Implementation of MetabR
A graphical user interface (GUI)based program, MetabR (Additional file 1), was written in the R opensource language (version 2.15). A screenshot of the GUI is shown in Figure 1. The GUI was built using the “gWidgets” package [24]. As described in the user guide (Additional file 2), the GUI is used to select which variables to define in a normalization model as fixed and randomeffect variables and to tailor statistical analysis to the researcher’s needs. As a threshold to screen for metabolites that differ significantly in abundance between treatment groups, the researcher may choose pvalue, qvalue, mean foldchange, or a combination of pvalue or qvalue and mean foldchange, as well as the specific values of these thresholds.
Additional file 1. MetabR. MetabR program file.
Format: R Size: 103KB Download file
Additional file 2. User Guide. MetabR user guide.
Format: PDF Size: 423KB Download file
This file can be viewed with: Adobe Acrobat Reader
Figure 1. Screenshot of the MetabR GUI.
In this program, either a fixed linear model (function “lm” in the “stats” package) or a linear mixed model (function “lmer” in the “lme4” package) [25] is constructed that includes the normalization variables selected by the user, and one Group treatment factor, for example
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.
To check normality and equal variance assumptions made by linear models, R functions “shapiro.test” in the “stats” package (“stats” and any other packages not referenced are part of R [26]) and “levene.test” in the “lawstat” package [27] are used, respectively. In addition, residual error plots are produced, and normalized data are exported for possible secondary use by the researcher. Tukey’s Honest Significant Difference (HSD; function “TukeyHSD” in the “stats” package) method is used to test for treatment group mean differences in the normalized data based on the Studentized range statistic. Qvalues [28] are calculated from the list of Tukey HSD pvalues for each treatment group comparison using the “qvalue” function in the “qvalue” package [28]. If any treatment group mean is significantly different from any other, a group mean plot with confidence intervals is constructed for the metabolite. Differences among treatment group means are represented by letter groupings generated by code adapted from a SAS macro [29], with means that share any letter being statistically equal. Further, statistical results from “significant” metabolites are exported into a spreadsheet that can be directly uploaded to Pathway Projector [9], which uses the information to map the metabolites, colored dots representing the direction and size of mean foldchanges, and either p or qvalues, to biochemical pathways. The program generates a series of files listed in Table 1 and described in the user guide (Additional file 2).
Table 1. Output files produced by the MetabR program
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 [30]. The extracted metabolomes were then analyzed by liquid chromatography—tandem mass spectrometry (LCMS/MS) via a slightly modified version of the methods of Rabinowitz and coworkers [3032] 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, signaltonoise 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 17dayold 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 files 3 (“Chicken example data 1”) and 4 (“Chicken example data 2”), corresponding to metabolites detected in positive and negative ionization modes, respectively.
Additional file 3. Chicken_pos. Chicken example data 1 from positive ionization mode.
Format: CSV Size: 4KB Download file
Additional file 4. Chicken_neg. Chicken example data 2 from negative ionization mode.
Format: CSV Size: 5KB Download file
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 16weekold 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 files 5 (“Mouse example data 1”) and 6 (“Mouse example data 2”), corresponding to metabolites detected in positive and negative ionization modes, respectively.
Additional file 5. Mouse_pos. Mouse example data 1 from positive ionization mode.
Format: CSV Size: 3KB Download file
Additional file 6. Mouse_neg. Mouse example data 2 from negative ionization mode.
Format: CSV Size: 3KB Download file
Modeling confounding variables as fixed vs. randomeffect
In our chicken example, Group, Quantity, and IS were modeled as fixedeffect variables, while Day was modeled as a randomeffect variable. To illustrate the difference, if Day is defined as a fixedeffect 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 randomeffect 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 fixedeffect variable, but experimental results can now be correctly extrapolated to all possible days [11].
Results
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 fixedeffect regression variables, and Day (run day) as a randomeffect factor.
Summary information printed in the R console (not shown) includes 1) results from the ShapiroWilk test of normality; 2) results from Levene’s test of equality of variance; 3) pairwise mean foldchanges between all treatment groups for significant metabolites (also exported into a spreadsheet; see Table 2 for the data); and 4) pairwise Tukey HSD pvalues or qvalues between all treatment groups for significant metabolites (also exported into a spreadsheet; see Table 2 for the data). This printout showed that ShapiroWilk pvalues for all metabolites were greater than 0.05, indicating no violation of the assumption of normality. Levene’s test pvalues for Citraconate and Inosine were less than 0.05, indicating a possible violation of the linear model assumption of equality of variance. By using the diagnostic results from ShapiroWilk and Levene’s tests, researchers can identify when data are unacceptable for use with the linear model approach. Ion counts are sometimes modeled as Poisson distributed, so if normality concerns are still an issue after opting for a log transformation, the researcher may wish to pursue an alternative statistical approach.
Table 2. Chicken experiment foldchanges
Figure 2 contains the plot of residual error for each metabolite after data transformation and normalization in relation to the overall mean abundance for each metabolite across all samples (we used log base 2 transformation). This plot can be used to determine whether data transformation and normalization corrected for the typical relationship of increasing variance with increasing mean. In this example, variance is visually relatively consistent across groups, and somewhat greater for lowabundance metabolites.
Figure 2. Residual error plot for the chicken experiment. Legend  Linear model residuals are plotted in relation to overall mean metabolite level.
Figure 3 illustrates an example of the mean plots and 95% confidence interval bars created for metabolites with a statistically significant effect of treatment. OAcetylLserine levels were significantly lower in Fast samples compared to InsNeut. Mean separation letters indicate that Fast and InsNeut groups differed significantly from each other (p < 0.05 threshold chosen), but neither differed from Control. Foldchanges between treatment group means (not log transformed) are displayed below the letters. Foldchanges in the nth row correspond to comparisons with the group in the nth column, (i.e., the mean of OAcetylLserine was 3.707fold higher in InsNeut compared to Fast).
Figure 3. Group mean plots for O AcetylLserine in the chicken experiment. Legend  Treatment group metabolite means, 95% confidence intervals, mean foldchanges, and significant difference letters are combined to summarize results for each significant metabolite.
Figure 4 shows a boxandwhisker plot of Citrate vs. run day, an ANOVA confounding variable, before and after data normalization with MetabR. Figure 5 shows a scatter plot of Pyruvate vs. Quantity, a regression confounding variable, before and after data normalization with MetabR. These plots are produced automatically by MetabR for all metabolites and all confounding variables included in the input data, and they give visual verification that the effects of confounding variables on metabolite measurements were removed by normalization using the linear model approach.
Figure 4. Pre and postnormalization plots: metabolite vs. Day. Legend  Citrate is plotted before and after normalization, showing the effectiveness of the normalization model for removing confounding variation in the chicken experiment. Normalization removed the effect of different run days on the Citrate detection signal.
Figure 5. Pre and postnormalization plots: metabolite vs. tissue quantity. Legend  Normalization removed the correlation between the quantity of tissue analyzed and the Pyruvate detection signal in the chicken experiment.
Figure 6 shows a heat map and dendrogram of the normalized data, produced automatically by MetabR via the “heatmap.2” function from the “gplots” package [35]. A heat map is useful for visualizing overall differences in metabolic signatures, and the dendrogram gives visual evidence of whether the experimental conditions significantly influenced metabolic signatures. Each metabolite plotted is meancentered, helping to call attention to metabolites differing in abundance among samples. The chickens appear to cluster nonrandomly based on their overall metabolic signatures (Figure 6). InsNeut chickens cluster in the upper half of the dendrogram, completely separate from Fast chickens, suggesting that these two treatment groups have distinct metabolic signatures, while the metabolic signature of the Control chickens appears less distinct.
Figure 6. Heat map and dendrogram. Legend  The heat map was produced by the MetabR program using the chicken example data included in Additional files 3 and 4. The metabolites are in the columns and the chicken adipose samples are in the rows. Columns are meancentered, with relative abundance represented by color (blue, lower abundance; red, higher abundance), as indicated in the legend. InsNeut chickens cluster in the upper half of the dendrogram, completely separate from Fast chickens, suggesting that these two treatment groups have distinct metabolic signatures, while the metabolic signature of the Control chickens appears less distinct. Note: the LCMS/MS instrument method is unable to differentiate between the several isomeric dihexoses, and therefore they are measured as a group.
Table 2 contains all betweengroup mean foldchanges for the metabolites, with differences tested by Tukey’s HSD at the 5% significance level. We produced this table by combining the mean foldchanges and pvalues exported automatically by MetabR. As shown, the experiment had sufficient power to detect a foldchange 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 qvalues automatically, and the researcher may select pvalue, qvalue, mean foldchange, or a combination of either pvalue or qvalue and mean foldchange 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 (qvalue) increasingly necessary.
Mouse experimental results
MetabR was run on the mouse example data in Additional files 5 and 6, selecting the same parameters as the chicken experiment, except that “Strain” (C57BL/6J or DBA/2J) was additionally selected as a randomeffect variable in order to remove the effects of different mouse strains on metabolite measurements. Two metabolites had ShapiroWilk pvalues less than 0.05 and W statistics less than 0.90, indicating possible violations of normality. No metabolites were identified as having unequal variance among treatment groups. The residual plot (not shown) also showed no evidence of unequal variance, and it was visually apparent that variance was equal across all measurement levels, and thus the log base 2 transformation chosen for the analysis was effective. Foldchange results are given in Table 3.
Table 3. Mouse experiment foldchanges
Conclusions
The opensource statistical computing software R [26] provides a convenient environment for statistical analysis of metabolomic and other omic data. We developed a userfriendly R program that normalizes metabolomic data using linear mixedeffect 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 RForge project at http://metabr.rforge.rproject.org/ webcite. 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.rforge.rproject.org/ webcite
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 nonacademics: No restrictions
Availability of supporting data
The datasets supporting the results of this article are included within the article (and its additional files).
Abbreviations
ANOVA: Analysis of variance; BPA: Bisphenol A; CSV: Commaseparated values; GUI: Graphical user interface; HSD: Honest Significant Difference; IS: Internal standard; LCMS: liquid chromatography—mass spectrometry; LCMS/MS: Liquid chromatography—tandem mass spectrometry.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
BE wrote the program. BE, SRC, BHV, and JRG collaborated to outline the issues in data analysis and process all biological data. AMS guided implementation of the statistical analysis components of the program. JRG and BE tested the implementation of the program, and all authors contributed to writing the final manuscript draft. All authors read and approved the final manuscript.
Authors’ information
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
Acknowledgements
JRG and SRC were supported by funding from the National Science Foundation through an Ocean Sciences award (OCE1061352) 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.
References

Nicholson JK, Connelly J, Lindon JC, Holmes E: Metabonomics: a platform for studying drug toxicity and gene function.
Nat Rev Drug Discov 2002, 1:153161. PubMed Abstract  Publisher Full Text

Reaves ML, Rabinowitz JD: Metabolomics in systems microbiology.
Curr Opin Biotechnol 2011, 22:1725. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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 AsianIndian men.
Diabetologia 2010, 53:757767. PubMed Abstract  Publisher Full Text

Kwon YKI, Higgins MB, Rabinowitz JD: Antifolateinduced depletion of intracellular glycine and purines inhibits thymineless death in E. coli.
ACS Chem Biol 2010, 5:787795. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Xia J, Psychogios N, Young N, Wishart DS: MetaboAnalyst: a web server for metabolomic data analysis and interpretation.
Nucleic Acids Res 2009, 37:W652W660. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Creek DJ, Jankevics A, Burgess KEV, Breitling R, Barrett MP: IDEOM: an Excel interface for analysis of LC–MSbased metabolomics data.
Bioinformatics 2012, 28:10481049. PubMed Abstract  Publisher Full Text

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:779787. PubMed Abstract  Publisher Full Text

Melamud E, Vastag L, Rabinowitz JD: Metabolomic analysis and visualization engine for LC − MS data.
Anal Chem 2010, 82:98189826. PubMed Abstract  Publisher Full Text

Kono N, Arakawa K, Ogawa R, Kido N, Oshita K, Ikegami K, Tamaki S, Tomita M: Pathway Projector: webbased zoomable pathway browser using KEGG Atlas and Google Maps API.
PLoS One 2009, 4:e7710. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Boccard J, Veuthey JL, Rudaz S: Knowledge discovery in metabolomics: an overview of MS data handling.
J Sep Sci 2010, 33:290304. PubMed Abstract  Publisher Full Text

Oberg L, Mahoney DH: Linear mixed effects models. In Topics in Biostatistics. Edited by Ambrosius WT. Totowa, NJ: Humana Press; 2007:213234.

Wolfinger RD, Gibson G: Assessing gene significance from cDNA microarray expression data via mixed models.
J Comput Biol 2001, 8:625637. PubMed Abstract  Publisher Full Text

Yang YH, Dudoit S, Luu P, Speed TP: Normalization for cDNA microarry data.

Berger MPF, Passos VL, Tan FES, Winkens B: Optimal designs for one and twocolor microarrays using mixed models: a comparative evaluation of their efficiencies.
J Comput Biol 2009, 16:6783. PubMed Abstract  Publisher Full Text

Chu TM, Weir B, Weir , Wolfinger R: A systematic statistical linear modeling approach to oligonucleotide array experiments.
Math Biosci 2002, 176:3551. PubMed Abstract  Publisher Full Text

Demirkale CY, Nettleton D, Maiti T: Linear mixed model selection for false discovery rate control in microarray data analysis.
Biometrics 2010, 66:621629. PubMed Abstract  Publisher Full Text

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: . PubMed Abstract  Publisher Full Text

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. BioMed Central Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Urs S, Smith C, Campbell B, Saxton AM, Taylor J, Zhang B, Snoddy J, Jones Voy B, MoustaidMoussa N: Gene expression profiling in human preadipocytes and adipocytes by microarray analysis.
J Nutr 2004, 134:762770. PubMed Abstract  Publisher Full Text

Wernisch L, Kendall SL, Soneji S, Wietzorrek A, Parish T, Hinds J, Butcher PD, Stoker NG: Analysis of wholegenome microarray replicates using mixed models.
Bioinformatics 2003, 19:5361. PubMed Abstract  Publisher Full Text

Smyth GK, Speed T: Normalization of cDNA microarray data.
Methods 2003, 31:265273. PubMed Abstract  Publisher Full Text

Du P, Kibbe WA, Lin SM: lumi: a pipeline for processing Illumina microarray.
Bioinformatics 2008, 24:15471548. PubMed Abstract  Publisher Full Text

Bates D, Maechler M, Bolker B:
lme4: Linear mixedeffects models using S4 classes. 2011.
[ http://CRAN.Rproject.org/package=lme4 webcite]

R Development Core Team: R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2011.
[ http:// www.Rproject.org/ webcite]

Noguchi K, Hui WW, Gel YR, Gastwirth JL, Miao W:
lawstat: An R package for biostatistics, public policy, and law. 2009.
[ http://CRAN.Rproject.org/package=lawstat webcite]

Storey JD: A Direct approach to false discovery rates.
Journal of the Royal Statistical Society B 2002, 64:479498. Publisher Full Text

Saxton AM: A macro for converting mean separation output to letter groupings in Proc Mixed. Nashville: Proceedings, 23rd SAS Users Group International: 2225 March 1998; 1996:12431246.

Collier JJ, Burke SJ, Eisenhauer ME, Lu D, Sapp RC, Frydman CJ, Campagna SR: Pancreatic βcell death in response to proinflammatory cytokines Is distinct from genuine apoptosis.
PLoS One 2011, 6:e22485. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bajad SU, Lu W, Kimball EH, Yuan J, Peterson C, Rabinowitz JD: Separation and quantitation of water soluble cellular metabolites by hydrophilic interaction chromatographytandem mass spectrometry.
Journal of Chromatography A 2006, 1125:7688. PubMed Abstract  Publisher Full Text

Waters CM, Lu W, Rabinowitz JD, Bassler BL: Quorum sensing controls biofilm formation in Vibrio cholerae through modulation of cyclic diGMP levels and repression of vpsT.
J Bacteriol 2008, 190:25272536. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Dupont J, Tesseraud S, Derouet M, Collin A, Rideau N, Crochet S, Godet E, CailleauAudouin E, MetayerCoustard S, Duclos MJ, Gespach C, Porter TE, Cogburn LA, Simon J: Insulin immunoneutralization in chicken: effects on insulin signaling and gene expression in liver and muscle.
J Endocrinol 2008, 197:531542. PubMed Abstract  Publisher Full Text

Ji B, Ernest B, Gooding J, Das S, Saxton A, Simon J, Dupont J, MetayerCoustard 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. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

gplots: Various R programming tools for plotting data. 2012.
[ http://CRAN.Rproject.org/package=gplots webcite]