Briefings in Bioinformatics Advance Access originally published online on October 31, 2006
Briefings in Bioinformatics 2007 8(2):71-77; doi:10.1093/bib/bbl019
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Enrichment analysis in high-throughput genomics—accounting for dependency in the NULL
Corresponding author. David L. Gold, Department of Statistics, Texas A&M University 3134 TAMU, College Station TX 77843-3143, USA. E-mail: dlgold{at}tamu.edu
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
Translating the overwhelming amount of data generated in high-throughput genomics experiments into biologically meaningful evidence, which may for example point to a series of biomarkers or hint at a relevant pathway, is a matter of great interest in bioinformatics these days. Genes showing similar experimental profiles, it is hypothesized, share biological mechanisms that if understood could provide clues to the molecular processes leading to pathological events. It is the topic of further study to learn if or how a priori information about the known genes may serve to explain coexpression.
One popular method of knowledge discovery in high-throughput genomics experiments, enrichment analysis (EA), seeks to infer if an interesting collection of genes is ‘enriched’ for a Consortium particular set of a priori Gene Ontology Consortium (GO) classes. For the purposes of statistical testing, the conventional methods offered in EA software implicitly assume independence between the GO classes. Genes may be annotated for more than one biological classification, and therefore the resulting test statistics of enrichment between GO classes can be highly dependent if the overlapping gene sets are relatively large. There is a need to formally determine if conventional EA results are robust to the independence assumption.
We derive the exact null distribution for testing enrichment of GO classes by relaxing the independence assumption using well-known statistical theory. In applications with publicly available data sets, our test results are similar to the conventional approach which assumes independence. We argue that the independence assumption is not detrimental.
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
The latest high-throughput genomics experiments generate an exhaustive quantity of measurements on typically thousands of genes. The analyst's role is not completed by producing an interesting list of genes associated with a phenotype or outcome. The results must be interpreted in a manner that makes biological sense. Working with higher-level information of common regulatory relationships or shared pathways, i.e. biological ‘classes’ of genes, shifts the focus to biological events that are easier to interpret than thousands of gene-wise events, and may further lead to hypotheses concerning specific genes. Inference on biological ‘classes’ of genes in high-throughput experiments has been referred to as functional enrichment, pathway analysis [1] or gene set enrichment analysis [2] (GSEA), but we refer to it simply as enrichment analysis (EA).
The gene ontology consortium [3] (GO) provides the most extensive and up-to-date prior knowledge of gene products available for high-throughput genomic analysis. GO is a 17-member organization that curates gene annotations for a diverse collection of model organisms ranging from protozoa to Homo sapien. The GO public access database, http://www.geneontology.org, consists of a controlled vocabulary describing biological classifications of genes embedded within three directed acyclic graphs, and a mapping from the known genes to the nodes on the graphs. Three separate graphs account for the biological processes (what), molecular functions (how) and cellular components (where) of the gene products. GO's hierarchically structured vocabulary provides multiple levels of description of gene products. For example, the child nodes under transcription in the biological process graph include: positive regulation of transcription and negative regulation of transcription. There are many publicly available GO browsers, such as the Cancer Genome Anatomy Project (CGAP) browser at http://cgap.nci.nih.gov/Genes/GOBrowser, which allow users to navigate through the dense vocabulary rooted within GO's three directed acyclic graphs.
In practice, EA is performed on pruned versions of the GO hierarchical trees. In their original form, many terminal leaf nodes on the GO trees provide a level of description that is too detailed for a higher-level analysis, while being associated with just a few genes. New terminal leaf nodes may be defined by specifying a maximum ‘depth’ for each branch; assigning all genes mapping to nodes beneath the maximum depth to the newly defined terminal leaf nodes. For example, if interest is generally in proliferation, genes associated with children of the cell proliferation node on the Biological Process graph may be reassigned to the cell proliferation parent class. Criterion for choosing a maximum tree depth for EA largely remains subjective despite statistical considerations, such as the number of GO classes to be tested for enrichment and the number of genes in and between each GO class.
|
FISHER EXACT TEST |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
The conventional approaches to EA test each GO class for enrichment individually, and report the resulting test statistic P-values without any adjustment to account for possible dependencies between the classes. The most popular method of inference is a one-way Fisher Exact Test. To fix ideas, consider a microarray experiment that measures the expression of N genes in multiple samples under different experimental conditions. Suppose that the preliminary analysis of these data produce a collection, C, containing k genes that are ‘interesting’ in the sense that they are either differentially expressed between two groups of samples or are associated in a significant way with outcome. Now let b1 be a known class of genes, such as a specific GO class; the problem we are concerned with is whether the collection C of interesting genes is significantly ‘enriched’ for genes belonging to the known class b1. Assume that N1 of the N unique genes on the array belong to the class b1. Assume further that (denoted in lower case) n1 of these genes belong to the collection C and N1 – n1 belong to the complement of C (i.e. do not belong to C).
Under the assumption that the k genes in C were sampled, independently and identically distributed (iid), uniformly at random without replacement from the N genes on the array, the marginal distribution of gene counts from collection C that belong to gene class b1 is hypergeometric. In other words, the probability of observing n1 genes from class b1 in collection C by chance is given by the formula
|
| (1) |
EA seeks to infer if the observed count n1 is so large as to cast doubt that it originated from (1). In EA, the null hypothesis (Ho) and alternative hypothesis (Ha) are often stated as
- Ho: The genes in C were sampled iid uniformly at random without replacement from the N genes on the array
- Ha: The genes in C were sampled from gene class b1 in greater proportion than N1/N; i.e. C was ‘enriched’ with genes from the b1 gene class.
- Ha: The genes in C were sampled from gene class b1 in greater proportion than N1/N; i.e. C was ‘enriched’ with genes from the b1 gene class.
The P-value for testing Ho is calculated by summing the probabilities in (1) from y = n1 to min (k, N1), which translates to measuring the probability of observing a count greater than or equal to n1, assuming that the null distribution in (1) is the true count distribution. A small P-value is evidence against Ho in favor of Ha. This procedure is essentially the one that R. A. Fisher [4] proposed for testing one-way independence of counts in 2 x 2 contingency tables, i.e. the one-way Fisher Exact Test.
The P-value cutoff for determining significance depends on the desired family-wise error rate (FWER). The FWER is the rate of detecting at least one false positive. In multiple testing, which occurs when a researcher tests for the enrichment of many different gene classes bi in the same collection C, the number of false positive test results is expected to grow with the number of tests. Classically, the P-value from a single test is adjusted in order to control the FWER. The most conservative approach is to multiply the single-test P-value by the number M of gene classes being tested, yielding the ‘Bonferroni-corrected’ P-value 
= pM. In other words, to ensure that FWER < , we must require the individual P-values to satisfy P < /M. Benjamini and Hochberg [5] discuss more modern approaches to the statistical issues involved in multiple testing, with special attention to methods that control the false discovery rate (FDR) instead of the FWER.
Statistical power, defined as the probability of detecting a true departure from the null hypothesis, is an important consideration for deciding which GO classes to include in EA. Figure 1 of our supplementary website (http://bib.oxfordjournals.org/) shows the statistical power of detecting enrichment for a single GO class, as a function of the fold increase in the expected count of genes in C relative to the expected count under Ho, for different total gene class sizes ranging from N1 = 25 to 300. The power of detecting enrichment for a 2-fold increase in the expected count under Ho is only modest for classes of between N1 = 25 and 100 total genes, while reasonably strong for N1 = 300 genes. In practice EA is performed on GO classes of diverse sizes. The total number of genes in each class is an important consideration for choosing GO classes to test. Unless there is prior biological information to justify inclusion, we recommend eliminating classes for which there is little statistical power of detecting enrichment. We also recommend eliminating classes with substantially overlapping gene sets. Including redundant classes can increase the false positive count, while yielding little additional information. We strongly advise working with unique genes rather than probes. The results from EA should not be arbitrarily based on the number of replicate or variant probes that target each gene, a product of the array design, but rather based on conclusions concerning which unique genes are biologically interesting.
Curtis et al. [1] provide an extensive list of available software, such as GO Miner and GOsurfer (accessible through dChip software), for performing EA using one-way Fisher Exact Tests. Many of these applications allow one to view a list of uploaded genes displayed on a GO directed acyclic graph under their respective nodes, with color codes to indicate up- or down-regulation. In our experience, these software packages report the P-values as though only a single statistical test were being performed, with no adjustment for multiple testing. Some software packages, such as GOsurfer implemented in DNA-Chip Analyzer Software (dCHIP), count probes rather than genes. Moreover, the underlying model does not account for dependence between gene classes that share substantial numbers of genes in common.
|
THE NULL DISTRIBUTION WITH DEPENDENCE |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
There is no restriction on the number of GO classes that a single unique gene may map to on the GO tree. In fact, many gene products are believed to perform multiple roles within a host of cellular events. GO classes may share genes, and the overlapping genes between GO classes can lead to positive correlation between the GO class counts in compelling collections of genes, even under the null hypothesis. Positive correlation between GO class counts may occur for other reasons, e.g. classes of genes that are related biologically, but typically we are not prepared to build such information into the null hypothesis. In this section, we derive the exact null distribution of gene counts for overlapping GO classes using well-known statistical theory, and offer a simulation approach for generating the null distribution as well as a normal approximation.
Let N be the number of unique genes on the array and let k be the number of genes in some biologically interesting collection C. Without loss of generality, suppose each gene is annotated with one or more of two GO classes. So, there are four possible annotations: b1 only, b2 only, both (b1 intersect b2) or neither (b1 union b2)c. Let N1 be the number of genes (out of the total number N) associated with b1, N2 with b2, and N3 with (b1 intersect b2). We denote the respective counts of genes, out of the k genes in the interesting collection C, using lower case: (n1, n2, n3). From an iid random sample without replacement of k genes from N, we are formally interested in the sampling distribution of y1 = n1 + n3 and y2 = n2 + n3. The joint distribution of (y1, y2) is
|
| (2) |
|
| (3) |
|
| (4) |
|
| (5) |
|
| (6) |
Notice that if (b1 intersect b2) is the empty set, then under (3) the joint distribution of (y1, y2) is expressed with N3 = 0 and n3 = 0, indicating negative covariance between the counts by (6) since k is finite. The negative covariance will tend to be small if k is small relative to N. In contrast, Fisher's Exact Test is applied without adjusting for overlap in the observed counts between GO classes. In other words, existing practice implicitly assumes that the joint distribution of (y1, y2) under the null distribution is
|
| (7) |
Computing exact P-values under the null for (2) is computationally cumbersome. An empirical approximation can be simulated given k, Ni and N, selecting k genes at random from N, and recording the counts of 
, for i = 1 to B classes and r = 1 to R replications. The empirical P-value is computed as the percent of simulated counts as or more extreme than those observed. Obtaining accurate empirical P-values according to the null joint distribution may also be computationally challenging, as controlling family-wise error rates requires setting exceedingly small P-value cutoffs with hundreds or even thousands of tests. We do not advocate simulation in practice, since the accuracy that EA requires in order to measure the frequency of extremely rare events would entail an impractical number of replications.
Alternatively, the null joint distribution given by formulas (2) and (3) may be approximated as multivariate normal (MVN), following from the multivariate extension of the Central Limit Theorem for proportions [4]. Under the null hypothesis, as the number of genes N tends to infinity, k/N 
(0, 1) and Ni/N Pi (0, 1) for i = 1, ... B, the joint distribution of the proportions of genes from C observed in each of the GO classes, pi = ni/k, tends asymptotically to a multivariate normal distribution with mean vector P = (P1, P2, ... PB,)T and covariance matrix 
,
|
| (8) |
may be acquired in closed form from formulas (5) and (6). The normal approximation depends on large k, Ni and N. Generally N is very large in array experiments. Certain combinations of k and Ni may be insufficient for the approximation in (9) to achieve a desired level of precision to its target standard normal distribution. Simulating from the null joint distribution can be helpful to determine minimum isotropic combinations of k and Ni required for the normal approximation to achieve a desired total error rate (sum of false positive and negative error rates) consistent with the target distribution.
The approximation in (8) provides a convenient theoretical shortcut to sidestep the tedious calculations required by (2) and (3), though in its present form lacks the potential to make individual inferences for each GO class. Now there exists a matrix D, such that the GO class proportions in (8) may be linearly transformed,
|
|
|
APPLICATION |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
We compared the results of our multivariate approximation, accounting for dependence, with Fisher's Exact test results, which do not formally adjust for dependence, using the candidate genes reported in three microarray studies comparing gene expression: (1) by breast cancer metastasis in two ER+/– groups [8], by type of treatment in human leukemia cells [9], and between multiple myeloma versus normal bone marrow tissue [10]. Each study reported a set of candidate genes believed to be changing with the experimental factors. Simulation results, not shown, with GO classes of different total sizes from the array experiments showed that the normal approximation to the null joint distribution after transformation in (9) achieved reasonable proximity to the standard normal distribution for GO classes of at least 300 genes. We tested GO classes with at least 300 genes in each study for comparison.
The total gene counts, along with the counts and analytic correlations between all pairs, of GO classes that we tested are reported on our supplementary website for all three studies (http://bib.oxfordjournals.org/). The GO classes with the highest analytic correlations tend to be the classes with the most genes where one class is a complete subset of the other. Histogram distributions of the analytic correlation coefficients between all pairs of GO classes in the Breast cancer study are shown on our supplementary website: Figure 2(a) for classes with at least 300 genes, and Figure 2(b) for classes with at least 50 genes. Correlation tends to be overall weaker among relatively smaller GO classes.
Tables 2–4
list the P-values from both methods for the top 10 classes, ranked by MVN P-value. The number of observed genes in each GO class, ni for i = 1, ..., 10, is listed along with the expected count of genes under the null hypothesis, k Ni/N. The orders of the P-values, whether accounting for correlation or not, are similar although the p-values from the MVN approximation tend to be slightly smaller in the Breast cancer study. The exact covariance matrices used to re-weight the counts in (9) are noticeably sparse.
|
|
|
We compared simulated P-values with the P-values obtained by Fisher and the MVN approximation for GO classes with greater than 50 genes in the Breast cancer study. Scatter plots of the P-values from all three methods are available in Figures 3(a–c) on our supplementary website (http://bib.oxfordjournals.org/). The simulated P-values were formed by generating 100 000 simulated counts for each GO class, transforming the simulated counts as in (9) in order to make independent inferences for each GO class. The smallest P-values largely agree in all three methods. The simulated versus MVN P-values tend to fall on the identity line, with a departure between 0.40 and 0.60 corresponding to relatively small GO classes.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
Knowledge discovery in bioinformatics is advancing steadily through the sharing and implementation of ideas across scientific disciplines. The vast quantities of information that must be synthesized, interpreted and reported in high-throughput genomic experiments is motivating discovery in previously unrelated fields, posing fresh challenges for creative thinking and testing old ideas. Success is attained by the diverse community that makes up bioinformatics, questioning preconceived notions, and challenging current ways of thinking.
EA is now a widely accepted and popular method for knowledge discovery in high-throughput genomics experiments. Its simplicity allows it to be generalized across platforms, and applied to a variety of experiments. We posed the question of robustness, challenging the conventional assumption of independence between GO classes, an assumption that we knew was incorrect, though might be advantageous. The results from the public experimental data that we analyzed show rank agreement between the MVN approximation, accounting correlation between classes, and the Fisher's Exact Test. We found that while there are GO classes that share many genes, on the whole the overlap is modest leading to sparse covariance matrices. In high-dimensional settings, if the true covariance matrix is sparse, then accounting for correlation is not expected to lead to dramatically different conclusion.
The normal approximation provides a convenient and fast tool for comparing the test results which account for correlation with the Fisher's Exact Test. By construction, the method supplies an interpretation (through the covariance matrix) for the relative performance of the two tests. While calculating exact significance levels to account for correlation is computationally burdensome, the normal approximation is quick to compute, although it requires large cluster and GO class count sizes. Simulating the null joint distribution can require an impractical number of replications in order to obtain empirical P-values of a desired level of accuracy to control the FWER or FDR, also presenting computational challenges. Simulation studies showed reasonable asymptotics for the MVN approximation with GO classes of at least 300 genes, and relative agreement with simulated P-values. The results based on the MVN approximation are suitable to conclude that the EA results of the Fisher's Exact Test are robust to the independence assumption for studies like those that we examined.
Annotations and data platforms are evolving rapidly. We encourage the Bioinformatics community to remain skeptical of the implicit assumptions made in conventional software packages, especially applied to ever larger future data sets. The choices of platform design and tests of hypotheses are choices that we should continue to consider carefully together as a research community. Designing high-throughput genomics experiments and methodologies that will lead to knowledge discovery is a system wide effort, benefiting from the shared input of collaboration at all levels.
Key Points
|
|
|
Acknowledgments |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
The research described here was at least in part funded by NCI grant CA104620.
|
FOOTNOTES |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
David L. Gold is currently a PhD candidate in the Department of Statistics at Texas A&M University.
Kevin R. Coombes is the Section Chief of Bioinformatics at M.D. Anderson Cancer Center. He received his PhD in Mathematics from The University of Chicago, IL.
Jing Wang is a statistician in the section of Bioinformatics at M.D. Anderson Cancer Center. He received his PhD in Biophysics, from the University of Manitoba, Winnipeg.
Bani Mallick is a full professor in the Department of Statistics at Texas A&M. He received his PhD in Statistics from the University of Connecticut.
Submitted: March 14, 2006. Received (in revised form): May 25, 2006.
|
References |
|---|
|
TOPABSTRACTINTRODUCTIONFISHER EXACT TESTTHE NULL DISTRIBUTION WITH...APPLICATIONDISCUSSIONFOOTNOTESAcknowledgmentsReferences |
|---|
- Curtis RK, Oresic M, Vidal-Puig A. Pathways to the analysis of microarray data. Trends in Biotechnology (2005) 23:429–35.[CrossRef][Web of Science][Medline]
- Mootha VK, Lindgren CM, Eriksson KF, et al. PGC-1alpha-responsive genes involved inoxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet (2003) 34:267–73.[CrossRef][Web of Science][Medline]
- The Gene Ontology Consortium. Gene Ontology: tool for the unification of biology. Nat Genet (2000) 25:25–9.[CrossRef][Web of Science][Medline]
- Agresti A. Categorical Data Analysis. (2002) New York: John Wiley & Sons Inc.
- Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: a Practical and Powerful Approach to Multiple Testing. J R Stat Soc B (1995) 57:289–300.
- Johnson NL, Kotz S, Balakrishnan N. Discrete Multivariate Distributions. (1997) John Wiley and Sons Inc.
- Gold DL, Coombes KR, Wang J, Mallick B. Testing Gene Class Enrichment in High-throughput Genomics UT MD Anderson Cancer Center Department of Biostatistics & Applied Mathematics Technical Report. (2005) New York. http://www.mdanderson.org/pdf/biostats_utmdabtr_004_05.pdf.
- Wang Y, Klijn JG, Zhang Y, et al. ‘Gene-expression to predict distant metastasis of lymph-node-negative primary breast cancer. Lancet (2005) 365:671–9.[Web of Science][Medline]
- Cheok MH, Yang W, Pui CH, et al. Treatment-specific changes in gene expression discriminate in vivo drug response in human leukemia cells. Nature Genet (2003) 34:85–90.[CrossRef][Web of Science][Medline]
- Zhan F, Hardin J, Kordsmeier B, et al. Global gene expression profiling of multiple myeloma, monoclonal gammopathy of undetermined significance, and normal bone marrow plasma cells. Blood (2002) 99:1745–57.
[Abstract/Free Full Text]
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  D. W. Huang, B. T. Sherman, and R. A. Lempicki Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists Nucleic Acids Res., January 1, 2009; 37(1): 1 - 13. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||