How to decide which are the most pertinent overly-represented features during gene set enrichment analysis

Uniform representation of data: DAGs defining sets partially ordered by the inclusion relation

We will denote S the set of objects considered in the remaining of this paper. For example, S can be the set of proteins of an organism. We consider that a neighborhood is a set N (of sets of elements of S) partially ordered by the inclusion relation ≺. Partially ordered sets (posets) are generally represented by Hasse diagrams in which there is an edge from y to x if and only if y covers x (denoted x ≺ y). This means that x ≺ y and there is no other element z such that x ≺ z ≺ y (see [8] for more details).

In the following, we consider that a neighborhood is a Hasse diagram (the DAG of figure 1b) that defines a poset (N, ≺).

Target sets pertinence

Several methods and tools use a similarity or dissimilarity index to compare sets, we can cite amongst others: FunSpec [9], BlastSets [7], GOStat [10], EASE [11], PANDORA [12], aBandApart [13], goCluster [14], see [1] for a review. Even though, these methods are using various similarity indices to compare the query and target sets (hypergeometric, binomial, χ², Fisher's exact test, or percentages), they all have in common that they consider only counts of elements such as the sizes of the query and target sets, or the number of common and differing elements. Thus, when comparing two sets, the bigger the number of common elements and the smaller the number of differing elements, the more similar they are considered.

Formally, this corresponds to any similarity index F between a query set Q and a target set T, such that F(Q, T) increases with |T ∩ Q| and decreases with |T - Q|. Then, given such a similarity index for the comparison of a given query set to a neighborhood, it is not necessary to perform the comparisons with all the target sets in the neighborhood. We introduce the notion of pertinence of a target set for its comparison to a given query set, which allows to consider target sets that are likely to have good similarity values (elements in common with the query) and to ignore target sets that will give redundant results (not different enough from other target sets because of the set composition dependencies in the Hasse diagram representing the neighborhood).

Our main observation is that redundancy is caused by two target sets when one includes the other and when they have either the same common elements or the same differing elements with the query. For example in figure 3, the target sets T₁ and T₃ are both redundant with T₂ for the query set Q. T₁ is redundant because it includes T₂ and has the same common elements (which implies that it has more differing elements) so it is less similar and it does not bring more information than T₂ alone. Similarly, T₃ is redundant because it is included in T₂ and has the same differing elements (i.e. less common elements). The fact that one target set includes the other is important for the biological meaning of the results. Let us consider T₂ and T₅ of figure 3. In this case, one should be tempted to decide that only T₂ is pertinent (same common elements and fewer differing elements). Actually, T₅ is also pertinent because the differing elements do not include those of T₂ and thus T₅ may be associated to a pertinent nonredundant biological meaning.

Mathematically, the observation above is written simply as follows:

Definition. A target set T in a neighborhood N is pertinent for its comparison to a given query set Q if and only if:

T ∩ Q ≠ ∅ (1)

∄T' ∈ N such that T' ⊂ T and T' ∩ Q = T ∩ Q (2)

∄T' ∈ N such that T' ⊂ T and T' - Q = T - Q (3)

This mathematical definition suggests that one must test all possible T' to decide if a target set T is pertinent. However, it is easy to deduce that only the parent and child nodes of T in the Hasse diagram representing the neighborhood N must be checked. Let us suppose that such a T' exists for (2) (resp. (3)). Then, due to the inclusion relation between T and T', all the sets in N on the path linking T and T' also satisfy (2) (resp. (3)), and then especially a child node (resp. a parent node) of T.

As only the parent and child nodes of T need to be considered, the test of pertinence can be performed on the number of common and differing elements. This is because if these numbers are equal then we are in presence of the same elements (inclusion relation).

As a result, the mathematical definition can be simplified into the following 3 rules (illustrated in figure 3) that are more suitable for the design of an efficient algorithm:

Rule 1: |T ∩ Q| ≠ 0

Rule 2: ∄T' such that T' ≺ T and |T ∩ Q| = |T' ∩ Q|

Rule 3: ∄T' such that T ≺ T' and |T - Q| = |T' - Q|

Algorithm for the identification and the comparison of pertinent target sets in the explicit representation

The pertinence rules allow us to define an algorithm (given in Algorithm 1 of figure 4) for the identification of target sets that are pertinent for their comparison to a given query set. Its principle is to search the DAG of the neighborhood, starting from the leaves corresponding to query elements (Rule 1) and exploring their ancestors to identify nodes satisfying the pertinence definition. This corresponds to a multiple sources breadth-first search, in which the queue is initialized with the nodes corresponding to the query elements. Each time a node is processed it is tested for pertinence. The search can stop at nodes including the query (Rule 2) or when the target set size is too big to give a significant similarity value. In the latter case, a test on the target set size is performed if an upper bound max_target_size can be computed theoretically (which is the case for most of the similarity models). In this algorithm, we suppose that the sets corresponding to the nodes of the DAG are available (which is not always feasible). The worst-case time complexity of a breadth-first search is O(V + E) where V is the number of vertices of the DAG and E is the number of edges. To test the pertinence of a node, we need (i) to compute the number of common and differing elements of the sets corresponding to the nodes, and (ii) to compare these values to the parent and child nodes to test if neither Rule 2 nor Rule 3 are violated. The computation of the number of common and differing elements for a node can be done in O(|S|), the maximum size of a set. The test of pertinence done in pertinent(Q, T) necessitates an access to all the parent and child nodes, which adds up to O(2E) supplementary tests. Thus, the worst-case time complexity of Algorithm 1 is O(|S|V + 3E) = O(|S|V + E). The worst-case time complexity occurs when all the nodes except the root include some but not all of the query elements. Let us consider the average case, in which we expect the query set to be small compared to the total number of elements |S|. The target sets sharing elements with the query represent only a subgraph of the DAG (figure 5b and 5c), and, the pertinent target sets should have sizes that are commensurate with the query size, which implies that they are deep in the DAG (figure 5c). Then, the number of nodes processed is typically very small compared to V. Moreover, the average target set size is small compared to |S|. Thus, the added |S| factor to the complexity may be considered as a constant and be negligible in the average case.

Algorithm 1 assumes that the set compositions are available. In the following, we introduce compact representations of neighborhoods and algorithms efficiently working on such representations.

Generic compact representation of neighborhoods

It is generally inefficient to explicitly store the composition of all the sets of a neighborhood. A compact representation is needed. Such a representation should permit one both to identify and to generate pertinent target sets efficiently, in a way that avoids the generation of all the sets and the traversal of the entire graph. Indeed, the uniform representation of neighborhoods by DAGs is adequate for compactness: we can store only the DAG defining the poset and reconstruct sets corresponding to nodes on the fly. In this compact representation, leaf nodes (nodes without successors) correspond to singleton sets (one for each element of S) and all the other nodes correspond to sets that can be built by searching the labels of reachable leaf nodes. For efficiency reasons, internal nodes are labeled with the size of the sets they represent as we will explain later. Figure 6 illustrates the compact representation corresponding to the DAG of figure 1b.

Algorithm for the identification of pertinent target sets in the generic compact representation

Like in Algorithm 1, the principle is to start from the leaves representing query elements and traverse the DAG to search for pertinent target sets among the sets sharing elements with the query. The difficulty we have to solve is that the set compositions are not available. We thus (i) store the size of the set corresponding to a node, as illustrated in figure 6, (ii) order the nodes in the queue by their size and (iii) propagate the common elements during the search.

In order to test Rule 2, we need the number of common elements of the node processed and its child nodes. With the breadth-first order, the node 'mRNA metabolic process' of size 6 of figure 6 would have been processed before the node of size 4 'RNA splicing' (shorter path from the leaves), and thus the common element b would not have been propagated yet. The solution is to maintain the queue ordered by the set sizes to ensure that sets smaller than the node processed (this includes all its descendants) have already been processed. This way, all the common elements have been propagated and Rule 2 can be tested at the level of the node processed.

In order to test Rule 3, we need the number of differing elements of the node processed and its parent nodes. Unfortunately, this number is not available at this time for the parent nodes because all the common elements may not have been propagated to the parent nodes yet. For example, the node 'RNA processing' (size 5) in figure 6 is processed before the element g has been propagated to the node 'RNA metabolic process' of size 7. The solution is to consider the node processed as the parent node and test if Rule 3 is not violated for its child nodes that is, we test the pertinence of its child nodes.

As a result, the pertinence decision is divided in two steps:

Step 1: Rule 2 is tested at the level of the node processed.

Step 2: Rule 3 is tested at the level of the child nodes of the node processed.

The efficiency of the resulting algorithm (given in Algorithm 2 in figure 7), compared to Algorithm 1, is only affected by the extraction of the next element of the queue, which must be ordered by the set sizes. As we can have at most |S| different sizes of set, the worst-case time complexity is affected by a factor of log |S| by using an adequate data structure. As for the previous algorithm, the average target set size is expected to be small compared to |S|, so the log |S| factor may be considered as a constant and be negligible.

Specific representations of neighborhoods

The previous compact representation is general and can be used for any neighborhood. Nonetheless, we identified cases where further time and space optimizations can be envisaged. It arises when:

In the following, we present efficient algorithms for the identification of pertinent target sets in these specific representations.

Algorithm for the identification and the comparison of pertinent target sets in the tree compact representation

The main advantage to searching pertinent target sets in a tree is that for a given node, the child nodes define a non overlapping partition of the set their parent node represents, and thus, only the number of common and differing elements need to be propagated.

Compared to the previous algorithm, this one avoids (i) the merging of the common elements for each node and (ii) the extraction of the next element of the queue. The tree is composed of |S| leaves and at most |S| - 1 nodes, thus, the worst-case time complexity of this algorithm is O(|S|).

Algorithm for the identification and the comparison of pertinent target sets in implicit compact representations

An implicit representation requires us to provide a specific algorithm for each different implied DAG. However, this loss in genericity allows considerable time saving in the search for pertinent target sets and considerable space savings for storing the neighborhoods. We have chosen to present the sets of adjacent genes on the chromosome because it can be described very simply and briefly, leads to a straightforward algorithm, and was often encountered in our experiments.

For our illustration, we only need to store the genes in the order they appear on the chromosome. Thus, the space requirement is θ(|S|), instead of θ(|S|²) needed for the DAG representation.

To identify pertinent target sets, we only need to know the position on the chromosome of each of the query elements. Then, each pair of positions defines a lower and an upper bound of an interval that, in turns, defines a set. For such a set to be pertinent, the bounds of the interval must be such that the position just before (resp. after) the lower (resp. upper) bound must not be an element of the query since it would violates Rule 3. Rule 2 holds because the bounds correspond to query elements. The worst-case time complexity of the resulting algorithm is O(|Q|²), Q being the query set. Compared to Algorithm 2 working on the generic compact representation, this algorithm spares (i) the merging of common elements and (ii) the extraction of the next element of the queue.

Query sets

The query sets of proteins correspond to protein complexes of the yeast Saccharomyces cerevisiae referenced at the MIPS [17] (version 14112005 filtered against a list of validated open reading frames from the GDR Genolevures [18]). The motivation for this choice is that once the proteins involved in a complex are identified, the next step is often to search for a molecular function or a biological process with which to annotate the newly grouped set of proteins. Moreover, the yeast proteome is well annotated with 4 211, 4 936 and 5 451 annotated gene products (among about 6 000) respectively for the molecular function, biological process and cellular component branch of the Gene Ontology (source: [19]). We extracted 1062 query sets of proteins, one for each protein complex.

Target sets/neighborhoods

This construction implies that when a protein is annotated with a GO term, all GO terms on the paths from this term to the root (more general terms) are also annotating this protein.

Space requirements

Efficiency

Redundancy reduction

Frequent itemset mining

The problem of identifying pertinent target sets resembles the frequent or closed itemset mining problem (see [21] for details) in many aspects. Unfortunately, the methods developed for frequent itemset mining cannot be applied to our context. Indeed, these methods rely on the anti-monotonic property of the support function (minimum frequency of the itemsets in the database). In our situation, the pertinence test is not anti-monotonic: a non pertinent target set can have ancestors that are pertinent. As a result, the pertinence definition cannot be used in the same way to prune the search. Moreover, closed itemsets permit the generation of all frequent itemsets contrary to pertinent target sets (all their subsets are not necessarily pertinent and may also not be present in the neighborhood).

Dissimilarity index based over-representation methods and tools

The methods we describe in this article improve the global quality of the results found by using a statistical test to decide the over-representation of a given feature in a given set (reviewed in [1]). Depending on the test performed (hypergeometric, binomial, χ²,...) and the correction for multiple testing (bonferroni, false discovery rate,...), the set of over-represented features will vary in size but the top features (very significant p-values) will remain essentially the same. Among those features, we clearly showed mathematically and also with biological results that a lot of them are actually redundant and non-informative. As it is based on the same principles (dissimilarity index and adjustment for multiple testing), our method can only perform at least as good as those others. For example, we submitted the query set of proteins of the complex 440.30.10 to GOStats [10] and we obtained 24 significantly enriched GO terms in the biological processes branch among which 7 are not in table 4. The observed differences (additional hits) are due to different versions of the data and to the multiple testing adjustment method (more low similarity hits). The additional hits are related to the pertinent hits found in table 4 and do not bring much additional insight to our query set biological function.

An original approach exploiting the GO structure was proposed in [22]. Like us, they consider the GO as a partially ordered set and work on the DAG. Their method and ours diverge due to the dissimilarity index. To score target sets, they define a pseudo distance which can be stated roughly as the average distance between the genes of the query and the target GO term. While this approach is formal and applicable to ontologies in general, it suffers some significant limitations. First computationally, for each query set they need to score all the GO terms. Second statistically, because of the use of a distance, they are only able to rank the GO terms and cannot assess the significance of the results. And finally, they also encounter the problem of redundancy and pertinence of the results. They partially address it by finding in the top ranked GO terms, the ones that are not comparable (i.e. the one that should bring more non redundant information).

Information content based methods

An interesting approach has been proposed by [23] to take into account the GO hierarchy. The difficulty to address when dealing with the GO hierarchy is that the level of a GO term in the hierarchy does not reflect the degree of specificity of this term. As a result, the degree of specificity (GO level) at which to look for enrichment should not be specified in the query because it can yield to misleading results or missed discoveries. In [23], the authors propose an information theoretic approach that allows to specify the degree of specificity desired for the enriched features. This is done by splitting the graph of the Gene Ontology into subgraphs. The split is such that the resulting subgraphs (partition of the GO terms) contain comparable information content, i.e. they concern the same number of genes. To illustrate their approach, they analyze the 'MAP00190 oxidative phosphorylation' set of proteins corresponding to GenMAPP proteins involved in oxidative phosphorylation. For this analysis, the Gene Ontology was split into 6 partitions among which a clear enrichment in 'transport' is visualized. In contrast, a corresponding GO biological process levelwise analysis performed at depth 2 exhibits visual enrichments in 'cellular process' and 'physiological process' which is misleading.

With our method, it is not possible to specify a given degree of specificity as with their tool GOPaD [23]. However, similar results can be obtained by constructing other neighborhoods for the Gene Ontology that would correspond to different level of specificity or information content. An alternative solution can also be to use GO Slim (reduced version of the GO aimed at giving an overview of its content) instead of the whole Gene Ontology. More generally, by looking at all the levels of the GO hierarchy, our method successfully identifies pertinent target sets, which automatically selects the most relevant levels to look at. Besides, our method is more generic in the sense that it can be applied to any hierarchically defined sets. For example, it can be applied to the hierarchical clustering of gene expression profiles which results in a dendogram (figure 2a) or the gene localization on chromosomes (figure 2b). In those cases, an information content method is of no help because the degree of specificity needs to be specified a priori and this is typically not known.

Integration of multiple data sources

Another trend in the field is the search for enrichment in combination of features by the use of multiple data sources. The direct approach consists in intersecting target sets as proposed in [24] where target sets of genes with composite GO annotations are obtained. This allows to find enrichments that are significant for the composite annotation (e.g. 'cation transport' and 'ATPase activity') while not being enriched in the original annotations (i.e. 'cation transport' alone or 'ATPase activity' alone). A similar approach has been proposed in [25] where frequent co-annotations (keywords, GO terms, and KEGG pathways) are mined. The principle is to search for frequent itemsets in the features of the query genes (features co-occurring frequently), and then to look at the significance of the enrichment in the combined features. Alternatively, the converse approach consists in the addition of GO term relationships such as is-involved-in as proposed in [26]. The principle is to augment the Gene Ontology to connect terms from different branches (sub-ontologies) to reflect the fact that a molecular function is involved in a biological process which takes place in a cellular component.

Although, the enrichment of combination of features is not addressed in this paper, similar results can be obtained by manipulating the neighborhoods. For example, it is possible to combine the GO biological process and the GO molecular function neighborhoods by adding nodes corresponding to the set intersections (composite annotations) to the Hasse diagram representing the neighborhood. Similarly, the augmentation of neighborhoods such as the additional Gene Ontology layer proposed in [26] can be achieved by adding the corresponding edges between GO terms and by propagating the gene products through the newly created paths. This could prove useful as we have seen in the results obtained for complex 440.30.10 with the gene YGL128c that the GO annotations are sometimes missing in a particular branch whereas present in another one.

Numerical features

Numerical features can be very interesting to consider for feature enrichment. For example, it can be used to discover that some genes of a query set are surprisingly close to each other on a chromosome, or that all the molecular weights of the query proteins fall within a surprisingly small range. To our knowledge, our approach is the sole capable of searching for enrichments in numerical features such as the gene localization on chromosomes (see figure 2a and the results section on implicit compact representations). This might be because (i) it is inefficient and sometimes unfeasible to store and compare all the sets corresponding to adjacent genes and (ii) because the redundancy in the results (if not filtered for pertinence) makes them unexploitable.