An Integrated Statistical Approach to Compare Transcriptomics Data across Experiments: A Case Study on the Identification of Candidate Target Genes of the Transcription Factor PPARα

Abstract

An effective strategy to elucidate the signal transduction cascades activated by a transcription factor is to compare the transcriptional profiles of wild type and transcription factor knockout models. Many statistical tests have been proposed for analyzing gene expression data, but most tests are based on pair-wise comparisons. Since the analysis of microarrays involves the testing of multiple hypotheses within one study, it is generally accepted that one should control for false positives by the false discovery rate (FDR). However, it has been reported that this may be an inappropriate metric for comparing data across different experiments. Here we propose an approach that addresses the above mentioned problem by the simultaneous testing and integration of the three hypotheses (contrasts) using the cell means ANOVA model. These three contrasts test for the effect of a treatment in wild type, gene knockout, and globally over all experimental groups. We illustrate our approach on microarray experiments that focused on the identification of candidate target genes and biological processes governed by the fatty acid sensing transcription factor PPARα in liver. Compared to the often applied FDR based across experiment comparison, our approach identified a conservative but less noisy set of candidate genes with same sensitivity and specificity. However, our method had the advantage of properly adjusting for multiple testing while integrating data from two experiments, and was driven by biological inference. Taken together, in this study we present a simple, yet efficient strategy to compare differential expression of genes across experiments while controlling for multiple hypothesis testing.

Keywords

false discovery rate gene expression transcriptomics microarray transcription factor peroxisome proliferator-activated receptor alpha ANOVA

Introduction

Genome-wide transcriptional profiling, or transcriptomics, is extensively used to study how cells respond to certain stimuli or to diagnose and predict clinical outcomes.^1–4 Transcription factors (TFs) are the key effectors which control gene expression. From a variety of research fields, including nutrition sciences, there is a major interest in characterizing the genes and networks that are controlled by transcription factors. Advances in genome-wide expression profiling methodologies and the availability of model systems offered new, powerful tools to address this.^5–11

An effective strategy to elucidate the signal transduction cascades activated by transcription factors is through transcriptional profiling. Transcription profiling can be applied on gain- and loss-of-function TF mutants, and changes in global gene expression that are associated with the various phenotypes could then be used for a comprehensive understanding of TF function.^6,7,11–13 To this end, transcription factor target genes have to be efficiently and accurately identified from the transcriptomics data set. It is important to realize that from a biological perspective, TF target genes are only those genes that do significantly respond upon treatment with a potent agonist or gain of function, in wild type but not mutant (knockout) models. However, from a statistical inference point of view the identification of biological relevant target genes from such 2 × 2 factorial experiments is less straight-forward.

It is generally accepted that statistical testing is required to reliable identify differentially expressed genes (reviewed in eg, Allison et al.¹⁴). Moreover, since the statistical analysis of microarrays involves the testing of multiple hypotheses (genes) within one study, it is necessary to control for false positives. A frequently used metric to quantify the level of confidence any particular gene is differentially expressed, that takes into account multiple testing, is the false discovery rate (FDR).¹⁴ Therefore in many studies a cutoff based on the FDR rather than P-value is used to select significantly regulated genes within experiments, which subsequently are compared across experiments to identify transcription factor target genes. However, Higdon et al¹⁵ reported that the use of the FDR and its associated q-value may result in inconsistent and misleading interpretation of the comparisons across different experiments, especially when the effect sizes of the experiments vary dramatically, as for example is the case when comparing effects of potent agonists in wild type and TF knockout models.

Therefore, the purpose of the work described in the current paper is to present a strategy that optimally integrates and controls for multiple hypothesis testing using data obtained from two biological systems that respond completely different to a treatment. We outline our approach using one of our datasets on the mouse peroxisome proliferator-activated receptor alpha (PPARα).¹⁶ PPARα is a TF belonging to the nuclear receptor superfamily, and is activated by a variety of compounds, including dietary fatty acids and their derivatives as well as synthetic agonists.^17–19

Material and Methods

Experimental data

We illustrate our approach (Fig. 1) on one of our publicly available datasets (Gene Expression Omnibus (GEO) accession: GSE8295). This dataset was generated to identify PPARα target genes in mouse liver,¹⁶ and was also used by Higdon et al¹⁵ to illustrate the inappropriateness of using the FDR as cut-off metric when comparing two transcriptomics experiments with different effect sizes.

Figure 1.

Overview of our integrated strategy.

Briefly, pure bred wild type (129S1/SvImJ) and PPARα-null (129S4/SvJae-Pparatm1Gonz/J) mice²⁰ were fed chow or chow supplemented with 0.1% WY14643 (Chemsyn, Lenexa, KS) for 5 days (n = 4 mice per group). WY14643, ({4-Chloro-6-[(2,3-dimethylphenyl)amino]-2-pyrimidinyl}-sulfanyl)acetic acid (CAS: 50892-23-4), is a chemical that was developed by the pharmaceutical industry to lower serum cholesterol. It is not used in clinical applications, but it is rather used as prototype chemical to induce peroxisome proliferation. WY14643 is a highly specific and potent agonist for PPARα and is therefore often used in studies on this nuclear receptor.^21,22 On the sixth day, mice were anaesthetized and livers were excised. Total RNA was prepared using TRIzol reagent (Invitrogen, Carlsbad, CA) followed by purification using the RNeasy mini kit (Qiagen, Hilden, Germany). RNA integrity was checked by chip analysis (Agilent 2100 bioanalyzer, Agilent Technologies, Amsterdam, The Netherlands) according to the manufacturer's instructions. RNA was judged as suitable for array hybridization only if samples exhibited intact bands corresponding to the 18S and 28S ribosomal RNA subunits, and displayed no chromosomal peaks or RNA degradation products, and had a RNA integrity number (RIN) above 8.0). The Affymetrix GeneChip RNA One cycle Amplification Kit (Afffymetrix, Santa Clara, CA) was used to prepare labeled cRNA from 5 μg of total RNA, which subsequently was hybridized on Affymetrix Mouse Genome 430 2.0 plus arrays. The animal study was approved by the Local Committee for Care and Use of Laboratory Animals.

Cell means ANOVA model

The dataset on the identification of PPARα target genes in mouse liver has a 2 × 2 factorial design; that is factor ‘treatment’ has 2 levels (WY, Control), as has the factor ‘genotype’ (wild type, knockout). Analysis of variance (ANOVA) is commonly used for analyzing data from experiments with multiple categorical factors.^23,24 To appropriately identify candidate PPARα target genes, we propose to perform and integrate three comparisons using the cell means ANOVA model.²⁵ For every probeset the model was defined as follows:

Y_{ijk} = μ_{ij} + ε_{ijk}

where Y_ijk is the expression of a probeset at i^th treatment (1 for WY, 2 for Control) in j^th strain of genotype (1 for WT, 2 for KO) and k^th replication (n = 4), μ is the mean value of i^th treatment and j^th strain of each gene, and ∊_ijk is a random error term which follows normal distribution with mean = 0 and variance = σ².

Formally, the definition of a contrast C is expressed below, using the notation μ_j for the j^th treatment mean:

{C=c}_{1} μ_{1} + c_{2} μ_{2} + … + c_{j} μ_{j} + … + c_{k} μ_{k}

where,

c_{1} + … + c_{j} + … + c_{k} = \sum_{j=1}^{k} c_{j} = 0

As stated before, from a biological perspective, candidate PPARα target genes are only those genes that do significantly respond upon treatment with the potent PPARα agonist WY14643 in wild type but not in PPARα knockout mice. Therefore three different contrasts (comparisons) from this 2 × 2 factorial experiment were combined to infer the probesets that were significantly and PPARα-dependently regulated. The different contrasts tested were (Table 1):

Table 1.

The contrasts defining the different hypotheses.

μ_ij	Levels	Contrast 1 H: μ₁₁ – μ₂₁ = 0	Contrast 2 H₀: μ₁₂ – μ₂₂ = 0	Global contrast H₀:(μ₁₁ – μ₂₁) – (μ₁₂ – μ₂₂) = 0
μ₁₁	WY, WT	1	0	1
μ₁₂	WY, KO	0	1	–1
μ₂₁	Con, WT	–1	0	–1
μ₂₂	Con, KO	0	–1	1

Contrast 1: H₀: μ₁₁ – μ₂₁ = 0 versus H₁: μ₁₁ – μ₂₁ ≠ 0, returning all probesets regulated in the wild type mice by the agonist WY;

Contrast 2: H₀: μ₁₂ – μ₂₂ = 0 versus H₁: μ₁₂ – μ₂₂ ≠ 0, returning all probesets regulated in the PPARα knockout mice by the agonist WY; and

Global Contrast: H₀:(μ₁₁ – μ₂₁) – (μ₁₂ – μ₂₂ = 0 versus H₁: (μ₁₁ – μ₂₁) – (μ₁₂ – μ₂₂) ≠ 0, returning the overall differential expressed probesets in wild type versus knockout mice groups after treatment with WY compared to control.

The PPARα-dependently regulated probesets were then identified by extracting those probesets that were only significantly regulated in both Contrast 1 and Global Contrast, and subsequently corrected for multiple testing.

Implementation

All analyses were performed in R,²⁶ using packages from the Bioconductor project.²⁷ Probesets were redefined according to Dai et al.²⁸ In this study, probes were reorganized based on Entrez Gene database, build 36, version 2 (remapped CDF version 12). Our workflow was as follows (note that since we used a remappped chip definition file based on the Entrez Gene database, the terms probeset and gene are used interchangeably):

Expression estimates were obtained by GC-robust multiarray (GCRMA) normalization, using the empirical Bayes approach to adjust background.²⁹

For each of the three above-mentioned contrasts, differentially expressed probesets (genes) were identified using linear models, as implemented in limma.³⁰ For each contrast probesets were selected based on P < 0.05.

Probesets that were common only in Contrast 1 and the combined Global Contrast were identified. This set of probesets represented only transcription factor regulated genes, and was designated X.

Multiple testing was corrected for using a false discovery rate method,³¹ based on the Global Contrast considering the output of all probesets. Probesets in X that satisfied the criterion FDR < 5% were considered to be transcription factor target genes.

An schematic overview of our implementation is also given in the Figure 1, and the R-code and other required files are available as supplemental material.

Validation

To validate our integrated approach, obtained results (Fig. 2) were compared to results from the across experiment comparison (Fig. 3) using two sets of well-established PPARα target genes obtained from a recent review (Table 1 from Rakhshandehroo et al³²). These sets are available as supplemental material. The true positive rate (sensitivity) as function of the false positive rate (1-specificity) for different cutoff points was plotted for both the across experiment comparison and our integrated approach using the R-library ROCR.³³ The partial area under the ROC curve was calculated using P = 0.2 (thus 1-specificity = 0.2) as cutoff. This cutoff value was chosen because for the identification of transcription factor target genes a high specificity is required (>80%) before considering its sensitivity.³⁴ In addition, the biological features that were overrepresented in the lists of candidate PPARα target genes that were generated on the basis of both approaches were analyzed with the software tool Ontologizer,³⁵ applying the ‘parent-child-union’ (PCU) algorithm and using the biological process ontology of Gene Ontology.

Figure 2.

Application of the integrated approach on PPARα dataset GSE8295.

Figure 3.

The FDR based across experiment comparison of PPARα dataset GSE8295.

Results and Discussion

Identification of candidate PPARα target genes

The application of transcriptomics to compare the effects of specific agonists, such as WY14643, in wild type and PPARα knockout mice is a powerful approach to identify candidate PPARα target genes.^12,13 However, when comparing across different experiments the use of FDR cutoff values may result in inconsistent and misleading interpretation of the data.¹⁵ In this study we propose a simple yet effective strategy that avoids comparing probesets across experiments based on FDR values while still controlling for multiple testing. Testing three different hypotheses (contrasts) for each probeset allowed the robust identification of transcription factor target genes. Since only the interaction effects are of interest for identifying candidate target genes, the cell means ANOVA model was used to infer this 2 × 2 factorial design.

The number of probesets significantly regulated (P < 0.05) upon PPARα activation by WY14643 in wild type mice (= Contrast 1) equaled to 5458, whereas in PPARα –/– mice (Contrast 2) this number was only 1540 (Fig. 2).

Such a large difference was expected since the KO mice do not express any functional PPARα. The Global Contrast, incorporating expression information for all probesets in all groups, identified 4282 significantly regulated probesets (P < 0.05) (Fig. 2), representing genes that from an inferential perspective are differentially regulated by WY between the two mouse strains. However, these included genes that for example were only regulated in the KO mice, or were regulated in wild type and, although to a lesser extent, still in KO mice. To filter out these ‘biological irrelevant’ genes, only probesets that were common in Contrast 1 and Global Contrast were retained, resulting in a set of 3345 probesets, which was called set X. Thus, this set X contained only probesets that from a biological perspective fulfill the criterion of being candidate PPARα target genes. To correct for multiple testing, FDR values of the 3345 genes in X were calculated based on all 16392 genes in Global Contrast, since in this comparison statistical inference was simultaneously adjusted for both experiments in wild type and knockout mice. Finally, a robust set of PPARα target genes was obtained by selecting those 2432 probesets from set X that fulfilled the criterion FDR < 0.05 (Fig. 2). Of these, 1325 probesets were induced and 1107 probesets were suppressed.

For comparison, we also generated a list of candidate PPARα target genes that were generated on the basis of directly comparing the wild type and knockout experiment using a FDR cutoff (Fig. 3). Note that this frequently used approach is criticized¹⁵ and that it is in essence identical to the analysis strategy published and interpreted by Rakhshandehroo et al¹⁶, except that these authors also employed a fold change cutoff. Using a FDR cutoff of 0.05, we identified 4324 probesets that were regulated in the wild type experiment (Contrast 1), whereas 12 probesets were changed in the knockout experiment (Contrast 2). Of these 12 probesets, 10 were also regulated in the wild type experiment. Thus, the FDR based comparison of these two experiments identified 4314 probesets that should be considered PPARα target genes.

The number of FDR based selected probesets was about twice as large as the list of probesets obtained using our integrated approach (4314 versus 2432 probesets). Comparison of these two sets of candidate genes revealed that almost all (i.e. 99%) of the probesets obtained by our integrated approach were also identified when using a FDR cutoff (Fig. 4). This indicates that while Global Contrast is more conservative it will identify similar if not identical biological features (see also section on validation).

Figure 4.

Venn diagram of the identified candidate PPARα target genes obtained by our integrated approach or the FDR based across experiment comparison of PPARα dataset GSE8295.

It is important to realize that the results of statistical hypothesis testing are never free of error. Two types of error are distinguished: type I error, i.e. rejecting the null hypothesis when it is in fact true, and type II error, i.e. not rejecting the null hypothesis when in fact the alternative hypothesis is true. In other words, occurrence of the former leads to inclusion of false positives whereas the latter leads to inclusion of false negatives. Consequently, we cannot exclude that the set of 1911 probesets that were discarded by Global Contrast contained false negatives that otherwise would have been retained. However, especially within the context of genome-wide screening studies for candidate genes, we believe that limiting type I error is of primary concern, and that of type II error is of secondary importance. Thus, to err on the safe side we prefer to control for false positives rather than for false negatives. Moreover, the probesets that were discarded by Global Contrast were characterized by a relatively low effect size compared to the probesets that were still included. The mean of the absolute coefficients (log2 of the fold-change) of the excluded probesets was 0.36 (equaling to a mean fold change of 1.28), and was 0.87 (mean FC = 1.83) for the included probesets. Taken together, we showed that compared to the FDR based across experiment comparison our approach identified a conservative set of more robustly regulated candidate PPARα target genes. We believe this is advantageous because a clear overview of candidate genes and corresponding biological process normally is aimed for.

Validation

To compare the performance of our integrated approach with that of the FDR based across experiments comparison, we first performed a sensitivity versus specificity analysis. To this end two benchmark sets of well-established PPARα target genes were selected from a review that summarized the latest literature on this topic.³² We created two benchmark sets; one set containing only 32 genes, and another set containing 189 genes. The smaller benchmark set contained only genes that were demonstrated to be bona fide PPARα target genes in both human and mouse liver and that do contain a functional PPAR response element (PPRE) in the regulatory regions. The larger benchmark set contained all genes that were demonstrated to be PPARα-dependently regulated in mouse liver but for which no functional PPRE has yet been identified. We next plotted the true positive rate (sensitivity) as function of the false positive rate (1-specificity) for different cutoff points for both our integrated approach and the across experiment comparison (Fig. 5).

Figure 5.

Sensitivity versus specificity of our proposed method and the across experiment comparison.

Even though our approached identified a conservative list of candidate genes, we observed that it performed very similar to the across experiment comparison in identifying known PPARα target genes, which was also reflected by almost identical partial area under the ROC curve (pAUC; P = 0.2) for both methods. Values were 0.129 and 0.121, respectively for the across experiment comparison and our integrated approach when the smaller set of 32 PPARα target genes was used, whereas these numbers were 0.128 and 0.124, respectively, for the larger set of 189 putative PPARα target genes.

Next we detected and compared the biological features that were overrepresented in the lists of candidate PPARα target genes that were either generated by our approach or the across experiment comparison. Enriched biological processes were identified by overrepresentation analysis based on Gene Ontology (GO) categories, which is a generally accepted procedure to achieve this.¹⁴ To this end the software tool Ontologizer was used,³⁵ applying the ‘parent-child-union’ (PCU) algorithm. This algorithm takes the graph structure of GO into account, thereby reducing false-positive and biologically misleading results.³⁶ We ran the PCU algorithm on the biological process ontology of GO. Using the Benjamini-Hochberg correction for multiple testing and a cutoff of 0.05, 27 significantly enriched categories (out of 172 annotated categories) were returned in the list of 4314 putative PPARα target genes obtained by comparing the wild type and knockout experiment (Fig. 6A).

Figure 6.

Significantly enriched Gene Ontology categories found in the two lists of candidate PPARα target genes.

Similarly, using the same criteria 28 significantly enriched categories (out of 169 annotated categories) were scored in the list of 2432 genes generated by our integrated approach (Fig. 6B). Twenty-five identified enriched biological processes were identical in both sets of genes. As expected, many processes that were enriched have been functionally demonstrated to be controlled by PPARα, including cellular ketone metabolic process, lipid metabolic process, cellular amino acid and derivative metabolic process, peroxisome organization, and mitochondrion organization.^17,37,38 Thus, despite the drastically reduced number of candidate PPARα target genes identified by our approach, GO enrichment analysis demonstrated a very similar functional characterization of these genes, again demonstrating the validity of our strategy.

Conclusions

Taken together, in this study we present a simple, yet efficient strategy to compare genes across experiments that controls for multiple testing and is able to properly detect differentially expressed genes. Compared to the conventional used FDR based across experiment comparison, our approach is more conservative and less noisy. Our approach is in particular suitable to identify candidate target genes of a transcription factor or signaling route from functional genomics experiments, but can be applied to any genomics experiment in which the effects of a treatment are compared between two genotypes.

Funding

This work was partially supported by the Wageningen University program for Systems Biology to MM (IP/OP SysBio).

Author Contributions

Conceived and designed the experiments: MOU, MM, GH. Analysed the data: MOU, GH. Wrote the first draft of the manuscript: MOU. Made critical revisions and approved final version: MM, GH. All authors reviewed and approved of the final manuscript.

Disclosures and Ethics

As a requirement of publication author(s) have provided to the publisher signed confirmation of compliance with legal and ethical obligations including but not limited to the following: authorship and contributorship, conflicts of interest, privacy and confidentiality and (where applicable) protection of human and animal research subjects. The authors have read and confirmed their agreement with the ICMJE authorship and conflict of interest criteria. The authors have also confirmed that this article is unique and not under consideration or published in any other publication, and that they have permission from rights holders to reproduce any copyrighted material. Any disclosures are made in this section. The external blind peer reviewers report no conflicts of interest.

Footnotes

The authors disclose no potential conflicts of interest.

Supplementary Materials

Supplementary 1. R code to perform the described analysis.

Supplementary 2. Sample description file for use with limma.

Supplementary 3. Excel file containing the two benchmark gene sets.

References

Brown

P.O.

, Botstein

Exploring the new world of the genome with DNA microarrays.

Nat Genet. 1999; 21(suppl 1): 33–7.

Stuart

J.M.

, Segal

, Koller

, Kim

SK.

A gene-coexpression network for global discovery of conserved genetic modules.

Science. 2003; 302(5643): 249–55.

Segal

, Friedman

, Kaminski

, Regev

, Koller

From signatures to models: understanding cancer using microarrays.

Nat Genet. 2005; 37 Suppl: S38–45.

Quackenbush

Microarray analysis and tumor classification.

N Engl J Med. 2006; 354(23): 2463–72.

Muller

, Kersten

Nutrigenomics: goals and strategies.

Nat Rev Genet. 2003; 4(4): 315–22.

Kirmizis

, Farnham

PJ.

Genomic approaches that aid in the identification of transcription factor target genes.

Exp Biol Med. 2004; 229(8): 705–21.

Gregory

B.D.

, Belostotsky

DA.

Whole-genome microarrays: applications and technical issues.

Methods Mol Biol. 2009; 553: 39–56.

Babu

M.M.

, Lang

, Aravind

Methods to reconstruct and compare transcriptional regulatory networks.

Methods Mol Biol. 2009; 541: 163–80.

Brien

G.L.

, Bracken

AP.

Transcriptomics: unravelling the biology of transcription factors and chromatin remodelers during development and differentiation.

Semin Cell Dev Biol. 2009; 20(7): 835–41.

10.

Petricka

J.J.

, Benfey

PN.

Reconstructing regulatory network transitions.

Trends Cell Biol. 2011; 21(8): 442–51.

11.

Gorte

, Horstman

, Page

R.B.

Microarray-based identification of transcription factor target genes.

Methods Mol. Biol. 2011; 754: 119–41.

12.

Bunger

, Hooiveld

G.J.

, Kersten

, Muller

Exploration of PPAR functions by microarray technology—a paradigm for nutrigenomics.

Biochim Biophys Acta. 2007; 1771(8): 1046–64.

13.

Woods

C.G.

, Heuvel

J.P.

, Rusyn

Genomic profiling in nuclear receptor-mediated toxicity.

Toxicol. Pathol. 2007; 35(4): 474–94.

14.

Allison

D.B.

, Cui

, Page

G.P.

, Sabripour

Microarray data analysis: from disarray to consolidation and consensus.

Nature Reviews Genetics. 2006; 7(1): 55–65.

15.

Higdon

, van Belle

, Kolker

A note on the false discovery rate and inconsistent comparisons between experiments.

Bioinformatics. 2008; 24(10): 1225–8.

16.

Rakhshandehroo

, Sanderson

L.M.

, Matilainen

Comprehensive analysis of PPARalpha-dependent regulation of hepatic lipid metabolism by expression profiling.

PPAR Res. 2007; 2007: 26839.

17.

Desvergne

, Wahli

Peroxisome proliferator-activated receptors: nuclear control of metabolism.

Endocr Rev. 1999; 20(5): 649–88.

18.

Sampath

, Ntambi

JM.

Polyunsaturated fatty acid regulation of gene expression.

Nutr Rev. 2004; 62(9): 333–9.

19.

Germain

, Staels

, Dacquet

, Spedding

, Laudet

Overview of nomenclature of nuclear receptors.

Pharmacol Rev. 2006; 58(4): 685–704.

20.

Lee

S.S.

, Pineau

, Drago

Targeted disruption of the alpha isoform of the peroxisome proliferator-activated receptor gene in mice results in abolishment of the pleiotropic effects of peroxisome proliferators.

Mol Cell Biol. 1995; 15(6): 3012–22.

21.

Issemann

, Green

Activation of a member of the steroid hormone receptor superfamily by peroxisome proliferators.

Nature. 1990; 347(6294): 645–50.

22.

Willson

T.M.

, Brown

P.J.

, Sternbach

D.D.

, Henke

BR.

The PPARs: from orphan receptors to drug discovery.

J Med Chem. 2000; 43(4): 527–50.

23.

Kerr

M.K.

, Martin

, Churchill

GA.

Analysis of variance for gene expression microarray data.

J Comput Biol. 2000; 7(6): 819–37.

24.

Ayroles

J.F.

, Gibson

Analysis of variance of microarray data.

Methods Enzymol. 2006; 411: 214–33.

25.

Kutner

M.H.

, Nachtsheim

C.J.

, Neter

, Li

Applied Linear Statistical Models. New York: McGraw-Hill, 2005.

26.

Ihaka

Gentleman RR: a language for data analysis and graphics.

J C omput Graph Stat. 1996; 5: 299–314.

27.

Gentleman

R.C.

, Carey

V.J.

, Bates

D.M.

Bioconductor: open software development for computational biology and bioinformatics.

Genome Biol. 2004; 5(10): R80.

28.

Dai

, Wang

, Boyd

A.D.

Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data.

Nucleic Acids Res. 2005; 33(20): e175.

29.

, Irizarry

R.A.

, Gentlemen

, Martinez-Murillo

, Spencer

A model-based background adjustment for oligonucleotide expression arrays.

J Am Stat Assoc. 2004; 99: 909–17.

30.

Smyth

limma: Linear Models for Microarray Data. In: Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Edited by Gentleman R, et al. New York: Springer New York; 2005: 397–420.

31.

Benjamini

, Hochberg

Controlling the false discovery rate: a practical and powerful approach to multiple testing.

J Royal Stat Soc Ser B. 1995; 57: 289–300.

32.

Rakhshandehroo

, Knoch

, Muller

, Kersten

Peroxisome proliferator-activated receptor alpha target genes.

PPAR Res. 2010.

33.

Sing

, Sander

, Beerenwinkel

, Lengauer

ROCR: visualizing classifier performance in R.

Bioinformatics. 2005; 21(20): 3940–1.

34.

Huber

, Scholtens

, Hahne

, Heydebreck

Differential Expression. In: Bioconductor Case Studies. Edited by Hahne F, et al. New York: Springer New York; 2008: 89–102.

35.

Bauer

, Grossmann

, Vingron

, Robinson

PN.

Ontologizer 2.0—a multifunctional tool for GO term enrichment analysis and data exploration.

Bioinformatics. 2008; 24(14): 1650–1.

36.

Grossmann

, Bauer

, Robinson

P.N.

, Vingron

Improved detection of overrepresentation of Gene-Ontology annotations with parent child analysis.

Bioinformatics. 2007; 23(22): 3024–31.

37.

Kersten

, Mandard

, Escher

The peroxisome proliferator-activated receptor alpha regulates amino acid metabolism.

FASEB J. 2001; 15(11): 1971–8.

38.

Mandard

, Muller

, Kersten

Peroxisome proliferator-activated receptor alpha target genes.

Cell Mol Life Sci. 2004; 61(4): 393–416.