ProbCD: enrichment analysis accounting for categorization uncertainty

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 over-representation or enhancement) analysis. In spite of efforts to create probabilistic annotations, especially in the Gene Ontology context, or to deal with uncertainty in high throughput-based 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 open-source R-based 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 on-line interface was created to allow usage by non-programmers and is available at: . 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 high-throughput experimental techniques and (ii) probabilistic gene annotation.

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 Fisher-like 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 multiple-test 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][9][10][11][12][13][14][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 protein-protein 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 high-throughput 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 high-throughput experimental techniques or (ii) probabilistic gene annotation.

Implementation
ProbCD is an open-source 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 user-friendly web-based interface for the software, which is not limited to any particular organism. The on-line interface and the source-code 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.
To motivate the general probabilistic model, it is useful to examine an arbitrary 2 × 2 example in the deterministic scenario: where all x's are the counts of a regular contingency table over the gene sets G and H. In its matrix representation: which is the generalization of the Bernoulli Process to more than two possible outcomes. The notation Z ~Be(p 1 ,ʜ, p n ) represents the distribution: The random variable X is a matrix representation of a k × 2 contingency table: where·is the usual dot-product, 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 non-categorical data, we recall the pivotal works by where 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 : 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: 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 ( ) 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 over-represented (or equivalently, the gene list is enriched for t) depending on user-defined 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. The point of the following illustration is to show that even ontology terms annotated with modest probabilities can be considered to be over-represented if the list of genes obtained behave in a supportive pattern. Consider a hypothetical organism with 100 genes annotated in several GO terms, as described in the Additional Files. The genes gene 1 to gene 20 are deterministically annotated to the ontology term t = a. In other words, assume that it is well known that these 20 genes have some given functionality a. The experiment, for example from a hypothetical proteomics dataset, yielded a deterministic list of differentially expressed (DE) genes ranging from gene 1 to gene 10 . The contingency table for this problem is, therefore: 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 probability-based 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 deterministic-based 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.
The point of the following illustration is to show that the incorporation of probabilistic annotation information does not always translate to addition of terms into the enrichment result, as in the example above, but it can also mean the exclusion of non-relevant terms. Consider a hypothetical organism with 1100 genes. Let the genes gene 1 to gene 100 be grouped together in a cluster H after some genomic sequence analysis. Let the term a be annotated deterministically (Additional Files) yielding the contingency table: In this situation, H is clearly enriched for a within any meaningful significance cutoff. Now let the same annota-tion 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 non-significant 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 down-rank potentially spurious enrichment results.
The following illustration shows that the use of ordered categories (k > 2) can produce useful results when additional information, regarding the order, is added. Consider a hypothetical organism with 4000 genes. In a hypothetical network analysis, let the genes be categorized, for simplicity and without loss of generality, in a deterministic fashion in a natural order: hubs (H), regular nodes (N) and leaves (L). Let the term a be annotated deterministically (Additional Files) yielding the contingency table: If one cannot express the difference between hubs and regular nodes in the enrichment analysis, the contingency table is forced to be described as: The most connected nodes in the network are not enriched for a considering the consolidated table above using either the Fisher's Exact Test (p-value = 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 up-regulated differentially expressed, not differentially expressed, and down-regulated 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 cell-cycle.
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 non-periodic 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 (p-value = 0.065).
Suspecting that this non-intuitive result could be due to the probability threshold chosen to select periodic genes, illustrated in the Figure 1, one could repeat the same analysis above building the contingency table considering the cutoffs ‫ސ‬ (gene i is periodic) ≥ 50%, 95%, 99% or 99.99%. The result of this repeated analysis is also non-intuitive since the p-values are: 0.12, 1.0, 1.0 and 1.0 for 50%, 95%, 99% and 99.99% cutoffs, respectively, meaning that increasing the stringency to define a gene as periodic only decreases the significance of the enrichment for GO:0007090.
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/S-specific 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: p-value of 0.047, 0.024, 0.0058, 0.086 and 0.024 for 50%, 75%, 95%, 99% and 99.99%, respectively.
The above analysis process is repeated for all GO terms, with the results available as Additional Files and summarized in Figure 2. This figure suggests that there is a large variability in the possible final outcome of an enrichment analysis depending on the probability cutoff used to build the associated contingency table. This variability is avoided by ProbCD because it directly takes into account the uncertainty in the data instead of introducing a discretization step (Figure 1). Figure 1 Probability of being periodic. The blue curve represents the probability of a gene being periodic (Pr) according to the model of [22]. The genes are sorted by probability values (rank) on the horizontal axis to facilitate the visualization. The red curve is the deterministic approximation using a 70% probability cutoff to consider a gene as periodic: ‫ސ‬ (gene i is periodic) ≥ 0.70 ⇒ ‫ސ‬ (gene i is periodic) = 1 and ‫ސ‬ (gene i is periodic) < 0.70 ⇒ ‫ސ‬ (gene i is periodic) = 0. This approximation labels 15% of the genes as periodic. Figure 3 shows that ProbCD considers more terms (vertical axis in Figure 3) containing the word "cell cycle", likely associated to periodically expressed genes, as significant if compared to the usual enrichment analysis in a wide range of significance values (p in Figure 3). Although this is not a proof, since one cannot be certain about which "cell cycle"-marked terms should be enriched, this is a reasonable indication that one can, in fact, avoid the discretization step when building the enrichment problem using ProbCD and obtain meaningful results.

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 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 context-dependent 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 pre-processing [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.
Fraction of "cell-cycle" GO terms selected as a function of the p-value Figure 3 Fraction of "cell-cycle" GO terms selected as a function of the p-value. The curves show the fraction of GO terms containing the word "cell-cycle" in their definition that are considered significant as a function of the significance cutoff (p). The red curve is obtained with ProbCD and all others are obtained with one of the probability cutoffs: 50%, 70%, 95%, 99% or 99.99%.
Venn diagram of over-represented terms