ProbCD: enrichment analysis accounting for categorization uncertainty
 Ricardo ZN Vêncio^{1} and
 Ilya Shmulevich^{1}Email author
DOI: 10.1186/147121058383
© Vêncio and Shmulevich; licensee BioMed Central Ltd. 2007
Received: 11 July 2007
Accepted: 12 October 2007
Published: 12 October 2007
Abstract
Background
As in many other areas of science, systems biology makes extensive use of statistical association and significance estimates in contingency tables, a type of categorical data analysis known in this field as enrichment (also overrepresentation or enhancement) analysis. In spite of efforts to create probabilistic annotations, especially in the Gene Ontology context, or to deal with uncertainty in high throughputbased datasets, current enrichment methods largely ignore this probabilistic information since they are mainly based on variants of the Fisher Exact Test.
Results
We developed an opensource Rbased software to deal with probabilistic categorical data analysis, ProbCD, that does not require a static contingency table. The contingency table for the enrichment problem is built using the expectation of a Bernoulli Scheme stochastic process given the categorization probabilities. An online interface was created to allow usage by nonprogrammers and is available at: http://xerad.systemsbiology.net/ProbCD/.
Conclusion
We present an analysis framework and software tools to address the issue of uncertainty in categorical data analysis. In particular, concerning the enrichment analysis, ProbCD can accommodate: (i) the stochastic nature of the highthroughput experimental techniques and (ii) probabilistic gene annotation.
Background
The systemlevel approach to data analysis known as enrichment analysis (also known as overrepresentation or enhancement analysis) is now commonplace. Moreover, the number of available software tools to perform such analysis is large (see [1, 2] for comprehensive reviews). The preferred way to formalize the enrichment problem is by means of a contingency table, often 2 × 2.
Besides measuring the statistical significance of the null hypothesis that the rows and columns are independent, as yielded by Fisher's Exact Test [3] and Fisherlike methods [1, 2], it is also possible to measure statistical association between a table's rows and columns [4] (a detailed discussion on significance vs. association in the enrichment problem context can be found in [5]).
Most of the attention in the enrichment analysis problem has focused on issues such as the search for the best multipletest correction or the implementation of better userfriendly software interfaces to facilitate biologist's exploratory work [1]. However, one of the limitations that the available approaches still share is that they assume, explicitly or implicitly, that one is able to construct the contingency table exactly, without uncertainty in populating its cells. Some efforts to consider ranked lists of genes, ranked by their reliability, were proposed to ameliorate the aforementioned limitations [6], however they do not work on the categorical data framework and incorpore the probabilitic information in a heuristic fashion [7].
Recently, the computational biology community has been witnessing an increasing interest in probabilistic approaches to gene annotation, particularly in the Gene Ontology (GO) context, as a realization of the limitations imposed by the traditional deterministic and contextindependent gene annotation schemes [8–15]. These efforts are motivated by: the necessity to assess the error propagation in automatic gene annotation [9, 15]; desire to include different types of evidence sources such as proteinprotein interaction [8, 13] or phylogenomics [10, 12] and annotation extrapolation from model organisms to others [11, 14]. Meanwhile, the probabilistic nature of data obtained by highthroughput measurement techniques is well recognized and a number of attempts to model it were proposed over the past decade in various experimental contexts [16, 17]. However, these efforts are not integrally taken into account when usual enrichment analysis is performed.
We describe a computational solution that is able to deal with the uncertainty introduced in enrichment analysis due to: (i) the stochastic nature of the results obtained with such highthroughput experimental techniques or (ii) probabilistic gene annotation.
Implementation
ProbCD is an opensource software designed to perform probabilistic categorical data analysis. ProbCD is written in R [18] with a level of modularity that makes it suitable to be incorporated by existing development efforts of integrative tools [19]. To facilitate the usage by researchers with no knowledge of R, we implemented a userfriendly webbased interface for the software, which is not limited to any particular organism. The online interface and the sourcecode are available on the project's website [20].
The idea behind ProbCD's implementation is to formally represent the intuitive process of building a contingency table in a probabilistic manner. Informally speaking, each element to be placed in the contingency table is not considered to be indivisible, but instead is "shared", according to probabilistic rules, among the contingency table's cells in a manner that is conceptually similar to fuzzy membership. The theoretical and computational implementation aspects are described in detail below.
Without loss of generality, the following descriptions are applied considering one particular ontology term t that is associated with a set of genes, named simply as G_{ t }. It should be noted that G_{ t }is not restricted to the Gene Ontology categorization and can be any kind of classification or annotation.
The vector q contains a probabilistic annotation for all g of the organism's genes: q_{ j }= ℙ (gene_{ j }∈ G_{ t }) for j ∈ {1,⋯, g}. This probabilistic annotation is assumed to be given, typically obtained from some analysis process. The deterministic scenario corresponds simply to ℙ (gene_{ j }∈ G_{ t }) ∈ {0, 1}, and hence is a special case.
The matrix P contains a probabilistic description for all k possible outcomes of the property being studied. Therefore, P is a k × g matrix with elements P_{i, j}= ℙ (gene_{ j }∈ outcome_{ i }) for j ∈ {1,⋯, g} and i ∈ {1,⋯, k}. This probabilistic description of the data uncertainty is assumed to be given.
where 1_{{}} is the indicator function.
Inspired by this representation, it is easy to see that the "hard" indicator functions may be substituted by Bernoulli random variables in order to account for the categorization uncertainty. Since all sets are finite, the indicator functions can be represented as vectors in {0, 1}^{ g }and the sums over all genes as dot products. In a generic scenario, with given nondeterministic P and q, the contingency table represented by XP, q is a random matrix that is difficult to describe in closed form. It is also not compatible with the statistical formalism supporting Fisher's Exact Test or other wellknown Fisherlike approaches, as these are not applicable to random tables.
where·is the usual dotproduct, a_{ i }= (A_{i, 1},⋯, A_{ i, g }) is a rowvector of a 2 × g binary matrix A such that (A_{1, j}, A_{2, j})q_{ j }~Be(q_{ j }, 1  q_{ j }) and d_{ i }= (D_{i, 1},⋯, D_{ i, g }) is a rowvector of a k × g binary matrix D such that (D_{1, j},⋯, D_{ k, j })(P_{1, j},⋯, P_{ k, j }) ~Be(P_{1, j},⋯, P_{ k, j }).
It is very easy to extend this framework for completely generic k × m tables (m > 2), but this would be outside the scope of the ontology enrichment problem.
To measure statistical association between rows and columns in contingency tables, analogously to correlations for noncategorical data, we recall the pivotal works by L.A. Goodman and W.H. Kruskal [4]. Depending on the problem under consideration, an appropriate association measure function ρ can be chosen. ProbCD calculates the statistical association accounting for the stochastic nature of the table's categorization, reporting ρ = ρ ($\mathbb{E}$ [XP, q]), where $\mathbb{E}$ is the expectation operator. If the categorical data is represented by a regular 2 × 2 matrix, then Yule's Q can be used as the statistical association function ρ ≡ Q : ℝ^{4} → [1, 1]. If one is dealing with ordered contingency tables, then GoodmanKruskal's gamma, ρ ≡ γ : ℝ^{2k}→ [1, 1], can be used since it is the generalization of Yule's Q. Considering nonordered categories, there is no analogy with the usual correlations in [1, 1] and in this case, as suggested by [4], Cramer's T is used with ρ ≡ T : ℝ^{2k}→ [0, 1].
All the association measures implemented can be calculated for $\mathbb{E}$ [XP, q] ∈ ℝ^{2k}, while 2 × 2 Fisher's Exact Test pvalue cannot, since it is a function in ℕ^{4} → [0, 1]. Moreover, a pvalue is related to the significance only, containing no information about the actual association level.
The dichotomous case, which is the simplest one, gives a more intuitive illustration on how the association is calculated in practice for the particular implementation: $\mathbb{E}$ [X_{1,1}P, q] = E_{1, 1} = P_{1, 1} q_{1} + ⋯ + P_{1, g}q_{ g }, $\mathbb{E}$ [X_{2, 1}P, q] = E_{2, 1} = (1  P_{1, 1}) q_{1} + ⋯ + (1  P_{1, g}) q_{ g }, $\mathbb{E}$ [X_{1, 2}P, q] = E_{1, 2} = P_{1, 1}(1  q_{1}) + ⋯ + P_{1, g}(1  q_{ g }), $\mathbb{E}$ [X_{2, 2}P, q] = E_{2, 2} = (1  P_{1, 1}) (1  q_{1}) + ⋯ + (1  P_{1, g}) (1  q_{ g }) and ρ = (E_{1, 1} E_{2, 2}  E_{1, 2} E_{2, 1})/(E_{1, 1} E_{2, 2} + E_{1, 2}E_{2, 1}), which corresponds to Yule's Q.
To measure the statistical significance of the estimated association, ProbCD uses a randomization approach. The null distribution for the association measure, ρ*, is proposed to be estimated from several permutation rounds. In each round a gene j receives randomly its probabilities (${P}_{1,j}^{\ast},\cdots ,{P}_{k,j}^{\ast}$) from one of the g possible columns of P and an association value is calculated. The significance of the statistical association between rows and columns in the contingency table is calculated as p = ℙ (ρ* ≥ ρ). A term t is significantly overrepresented (or equivalently, the gene list is enriched for t) depending on userdefined thresholds for significance and/or association.
Results
The following examples illustrate the potential utility of considering probabilistic annotations and/or data uncertainty assessment in the enrichment analysis using ProbCD on artificial datasets and a published yeast dataset.
In this case, the DE gene list is clearly enriched for a within any meaningful significance cutoff. Consider now a second ontology term b obtained from a probabilitybased source with ℙ (gene_{ i }∈ G_{ b }) = 40%, i ∈ {1,⋯, 20}. A probability of only 40% generally would not be sufficient evidence to warrant the inclusion of those 20 genes in G_{ b }considering a usual deterministic framework and, therefore, would not be analyzed by deterministicbased methods, such as the Fisher's Exact Test. However, ProbCD is able to incorporate this information and yields: ρ = 0.87 and p < 10^{4} in 10000 permutation rounds, a significant enrichment for b. One can easily imagine, for example, genes that have a main function a but also have a different function b in, say, 40% of documented conditions.
In this situation, H is clearly enriched for a within any meaningful significance cutoff. Now let the same annotation incorporate some evidence levels by defining: ℙ (gene_{ i }∈ G_{ a }) = 99% for i ∈ {1,⋯, 10} and ℙ (gene_{ i }∈ G_{ a }) = 1% for i ∈ {11,⋯, 100}. Intuitively, this means that only 10 out of 100 genes clustered in H are, in fact, confidently annotated with the ontology term a. The incorporation of this information results in nonsignificant enrichment of H for a since: ρ = 0.0425 and p = 0.42 in 1000 permutation rounds. Therefore, it can be useful to incorporate uncertainty information into the enrichment analysis to also downrank potentially spurious enrichment results.
The most connected nodes in the network are not enriched for a considering the consolidated table above using either the Fisher's Exact Test (pvalue = 0.54) or ProbCD (ρ = 0, p = 0.48). However, using the original categorization order, ProbCD suggests a significant enrichment for a with ρ = 0.98 and p < 10^{4}. The conclusion that the property a must be related to gene connectivity seems subjectively reasonable considering the numbers in the first contingency table. The rationale used for the hypothetical network analysis could be useful in other scenarios where there is a natural order that can provide extra information such as: highly expressed, expressed, and not expressed or upregulated differentially expressed, not differentially expressed, and downregulated differentially expressed.
The next illustration demonstrates the impact of considering the uncertainty in lists of genes, rather than in the annotations, on the enrichment analysis. In this example, the aim is to find which GO terms, annotating the yeast Saccharomyces cerevisiae, are statistically associated with periodic expression levels, measured by microarray technology [22]. Andersson and colleagues [22] devised an elaborate Bayesian model which produces the probability that a gene is periodically expressed during the cellcycle. Since the final probability values are sufficient for our objectives in this work, we refer the interested reader to the original work by Andersson and colleagues [22] for more details. In this example, the annotation is considered to be deterministic and was downloaded from the GO project page (March 2007) [23].
To perform the usual enrichment analysis one needs to define a probability cutoff value in order to split the gene list in two: the periodic genes and the nonperiodic genes. Consider initially the reasonable cutoff ℙ (gene_{ i }is periodic) ≥ 70% and focus on a single GO term GO:0007090 (regulation of S phase of mitotic cell cycle), defined as "a cell cycle process that modulates the frequency, rate or extent of the progression through the S phase of mitotic cell cycle". Although this GO term is clearly associated with periodic gene expression, performing a usual enrichment analysis results in the conclusion that the periodic genes are not significantly enriched for GO:0007090 within usual significance cutoffs (pvalue = 0.065).
Using ProbCD, one can consider the actual probability of being periodic (blue curve in Figure 1) in the enrichment analysis instead of using the deterministic approximation (red curve in Figure 1). This results in a relatively high statistical association between periodicity and the term "regulation of S phase of mitotic cell cycle" (ρ = 0.78) with high significance (p = 0.009 in 1000 simulation rounds). Judging subjectively by the definition of GO:0007090, ProbCD returned a meaningful result.
Other similar cases can be easily identified. For example, the GO term GO:0000083 (G1/Sspecific transcription in mitotic cell cycle) exhibits erratic behavior depending on the chosen cutoff for the probability of being periodic: pvalue of 0.15, 0.10, 0.01, 0.096 and 1.0 for 50%, 75%, 95%, 99% and 99.99% cutoffs, respectively. The probability stringency used to build the contingency table and the subsequent significance test are not necessarily correlated. ProbCD yielded a significant (p = 0.006) moderate association (ρ = 0.48) for GO:0000083. Other examples include GO:0045787 (positive regulation of progression through cell cycle), defined as "any process that activates or increases the frequency, rate or extent of progression through the cell cycle", which would be called significant using the regular enrichment method only if the right probability cutoff ℙ (gene_{ i }is periodic) ≥ 95% is guessed initially: pvalue of 0.047, 0.024, 0.0058, 0.086 and 0.024 for 50%, 75%, 95%, 99% and 99.99%, respectively.
Discussion and Conclusion
The usual enrichment analysis is a particular case in this probabilistic framework and can be obtained by ProbCD ignoring the difference between evidence sources in gene annotation and defining fixed gene lists, which would correspond to the deterministic setting: q_{ j }= ℙ (gene_{ j }∈ G_{ t }) = 1 or 0 and P_{ i, j }= ℙ (gene_{ j }∈ outcome_{ i }) = 1 or 0.
Even if a probabilistic annotation is not readily available for a given organism, it could be interesting to perform enrichment analysis taking into account some form of weighting on available annotations according to their reliability. For a concrete example, the GO Consortium [24] provides annotations accompanied with evidence codes related to the kind/level of evidence available for a given GO annotation [25], such as IEA: Inferred from Electronic Annotation, IMP: Inferred from Mutant Phenotype, RCA: inferred from Reviewed Computational Analysis or IDA: Inferred from Direct Assay. It is known that some evidence sources are more reliable than others and this knowledge can be used, in a Bayesian sense, as subjective probabilities.
Once an annotation is considered in a probabilistic framework, it could reflect a dependence on the context. One can consider cases in which ℙ (gene_{ j }∈ G_{ t }disease) ≫ ℙ (gene_{ j }∈ G_{ t }), defining contextdependent gene annotations derived, for instance, from automatic literature mining [26].
Our intention is to complement existing approaches, rather than substitute them. Toward this aim, we built ProbCD to be as modular as possible in order to be incorporated into existent software or pipelines [19], composed of ontology preprocessing [27] or powerful visualization capabilities [28, 29].
It is important to note that ProbCD is also applicable to other categorical data analysis contexts in which the construction of contingency tables is subject to uncertainty, a recurrent theme in science.
Availability and requirements

Project Name: ProbCD

Project Home Page: http://xerad.systemsbiology.net/ProbCD

Operating Systems: platform independent

Programming Languages: R

License: GNU Lesser General Public License 3.0
Declarations
Acknowledgements
We thank Drs. John Boyle, Vesteinn Thorsson, Nathan Price and other ISB colleagues for insightful discussions. We also thank Dr. Carlos A. de B. Pereira (USP) for introducing us to the works of L.A. Goodman and W.H. Kruskal. This work is partially supported by NIH grants U54AI54253, U19AI057266 and P50GMO76547.
Authors’ Affiliations
References
 Dopazo J: Functional Interpretation of Microarray Experiments. OMICS: A Journal of Integrative Biology 2006., 10(3):Google Scholar
 Rivals I, Personnaz L, Taing L, Potier M: Enrichment or depletion of a GO category within a class of genes: which test? Bioinformatics 2007, 23(4):401–407.View ArticlePubMedGoogle Scholar
 Fisher R: On the Interpretation of χ^{2}from Contingency Tables, and the Calculation of P. Journal of the Royal Statistical Society 1922, 85: 87–94.View ArticleGoogle Scholar
 Goodman L, Kruskal W: Measures of Association for Cross Classifications. Journal of the American Statistical Association 1954, 49(268):732–764.Google Scholar
 Vencio R, Koide T, Gomes S, Pereira C: BayGO: Bayesian analysis of ontology term enrichment in microarray data. BMC Bioinformatics 2006, 7: 86.PubMed CentralView ArticlePubMedGoogle Scholar
 Jiang Z, Gentleman R: Extensions to gene set enrichment. Bioinformatics 2007, 23(3):306.View ArticlePubMedGoogle Scholar
 Goeman J, Buhlmann P: Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 2007, 23(8):980.View ArticlePubMedGoogle Scholar
 Joshi T, Chen Y, Becker J, Alexandrov N, Xu D: GenomeScale Gene Function Prediction Using Multiple Sources of HighThroughput Data in Yeast Saccharomyces cerevisiae. Omics A Journal of Integrative Biology 2004, 8(4):322–333.View ArticlePubMedGoogle Scholar
 Levy E, Ouzounis C, Gilks W, Audit B: Probabilistic annotation of protein sequences based on functional classifications. BMC Bioinformatics 2005, 6: 302.PubMed CentralView ArticlePubMedGoogle Scholar
 Engelhardt B, Jordan M, Muratore K, Brenner S: Protein molecular function prediction by Bayesian phylogenomics. PLoS Comput Biol 2005., 1(5):Google Scholar
 Martin D, Berriman M, Barton G: GOtcha: a new method for prediction of protein function assessed by the annotation of seven genomes. BMC Bioinformatics 2004, 5: 178.PubMed CentralView ArticlePubMedGoogle Scholar
 Engelhardt B, Jordan M, Brenner S: A graphical model for predicting protein molecular function. Proceedings of the 23rd international conference on Machine learning 2006, 297–304.Google Scholar
 Carroll S, Pavlovic V: Protein classification using probabilistic chain graphs and the Gene Ontology structure. Bioinformatics 2006, 22(15):1871.View ArticlePubMedGoogle Scholar
 Vinayagam A, del Val C, Schubert F, Eils R, Glatting K, Suhai S, König R: GOPET: A tool for automated predictions of Gene Ontology terms. BMC Bioinformatics 2006, 7: 161.PubMed CentralView ArticlePubMedGoogle Scholar
 Jones C, Brown A, Baumann U: Estimating the annotation error rate of curated GO database sequence annotations. BMC Bioinformatics 2007, 8: 170.PubMed CentralView ArticlePubMedGoogle Scholar
 Zhang W, Shmulevich I: Computational and Statistical Approaches to Genomics. 2nd edition. New York, NY, USA: Springer; 2006.View ArticleGoogle Scholar
 Zhang W, Shmulevich I, Astola J: Microarray Quality Control. WileyLiss; 2004.View ArticleGoogle Scholar
 The R Project for Statistical Computing[http://www.rproject.org]
 Shannon P, Reiss D, Bonneau R, Baliga N: Gaggle: An opensource software system for integrating bioinformatics software and data sources. BMC Bioinformatics 2006, 7: 176.PubMed CentralView ArticlePubMedGoogle Scholar
 ProbCD Home Page[http://xerad.systemsbiology.net/ProbCD]
 Bernoulli scheme – Wikipedia, The Free Encyclopediahttps://en.wikipedia.org/w/index.php?title=Bernoulli%20scheme&o%25ldid=64557593
 Andersson C, Isaksson A, Gustafsson M: Bayesian detection of periodic mRNA time profiles without use of training examples. BMC Bioinformatics 2006, 7: 63.PubMed CentralView ArticlePubMedGoogle Scholar
 Gene Ontology Current Annotations[http://www.geneontology.org/GO.current.annotations.shtml]
 The Gene Ontology Consortium[http://www.geneontology.org]
 Guide to GO Evidence Codes[http://www.geneontology.org/GO.evidence.shtml]
 Aubry M, Monnier A, Chicault C, de Tayrac M, Galibert M, Burgun A, Mosser J: Combining evidence, biomedical literature and statistical dependence: new insights for functional annotation of gene sets. BMC Bioinformatics 2006, 7: 241.PubMed CentralView ArticlePubMedGoogle Scholar
 Lewin A, Grieve I: Grouping Gene Ontology terms to improve the assessment of gene set enrichment in microarray data. BMC Bioinformatics 2006, 7: 426.PubMed CentralView ArticlePubMedGoogle Scholar
 Maere S, Heymans K, Kuiper M: BiNGO: a Cytoscape plugin to assess overrepresentation of Gene Ontology categories in Biological Networks. Bioinformatics 2005, 21(16):3448–3449.View ArticlePubMedGoogle Scholar
 Sealfon R, Hibbs M, Huttenhower C, Myers C, Troyanskaya O: GOLEM: an interactive graphbased geneontology navigation and analysis tool. BMC Bioinformatics 2006, 7: 443.PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
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.