Enumerating consistent subgraphs of directed acyclic graphs: an insight into biomedical ontologies
Yisu Peng Yuxiang Jiang Predrag Radivojac
cs.DS cs.DM q-bio.BM









Department of Computer Science,
 Indiana University, Bloomington, Indiana, U.S.A.

Enumerating consistent subgraphs of directed acyclic graphs: an insight into biomedical ontologies
    Predrag Radivojac[Corresponding author, email: ]
    December 30, 2023
==================================================================================================




Modern problems of concept annotation associate an object of interest (gene, individual, text document) with a set of interrelated textual descriptors (functions, diseases, topics), often organized in concept hierarchies or ontologies. Most ontologies can be seen as directed acyclic graphs, where nodes represent concepts and edges represent relational ties between these concepts. Given an ontology graph, each object can only be annotated by a consistent subgraph; that is, a subgraph such that if an object is annotated by a particular concept, it must also be annotated by all other concepts that generalize it. Ontologies therefore provide a compact representation of a large space of possible consistent subgraphs; however, until now we have not been aware of a practical algorithm that can enumerate such annotation spaces for a given ontology. In this work we propose an algorithm for enumerating consistent subgraphs of directed acyclic graphs. The algorithm recursively partitions the graph into strictly smaller graphs until the resulting graph becomes a rooted tree (forest), for which a linear-time solution is computed. It then combines the tallies from graphs created in the recursion to obtain the final count. We prove the correctness of this algorithm and then apply it to characterize four major biomedical ontologies. We believe this work provides valuable insights into concept annotation spaces and predictability of ontological annotation.




§ INTRODUCTION


Ontologies have become a common means of concept annotation in computational biology and related fields <cit.>. A protein's molecular function <cit.>, an effect of a genetic variant <cit.>, or a patient's diagnosis <cit.> are typical examples in which biomedical entities such as macromolecules, mutations, or individuals are associated with sets of mutually dependent descriptors. The dependencies between these descriptors are often hierarchical, leading to the use of directed acyclic graphs as concept space representations.

A directed acyclic graph is a pair (V,E), where V is a set of vertices (nodes) and E is a set of directed edges (links) between vertices such that no cycles can be formed. Each vertex in the graph is associated with a unique concept (term, description) and each edge is associated with a particular type of relational tie. For example, when annotating proteins as biomedical entities using the Gene Ontology graph <cit.>, the terms “nucleic acid binding” and “DNA binding” are linked by edges of the type  asserting that DNA binding is a more specific form of nucleic acid binding. Other types of relational ties include , , and so on.

A typical biomedical entity is associated with a set of terms determined through experiment such as a molecular assay or a diagnostic procedure. A protein, for example, may be assigned terms “DNA binding” and “RNA binding”, neither of which is a generalization of the other. To avoid annotation inconsistencies, this protein must also be annotated by the terms such as “nucleic acid binding” and all others that generalize either of the experimentally determined terms. More broadly, this implies that a biomedical object can only be annotated by a set of terms that respect the hierarchy – a consistent subgraph of the ontology. Unfortunately, (manual) experimental annotation is resource-demanding and often incomplete <cit.>, giving rise to an entire field of computational prediction <cit.>. 

The development of computational prediction methods presents its own challenges. Although it can be performed by building a separate binary classifier for each concept in the ontology, this approach is currently competitive only for specialized ranking tasks; e.g., disease-gene prioritization <cit.>, since it does not exploit relationships between the terms. On the other hand, a more complete characterization is via learning structured outputs <cit.> in which a method takes an object (e.g., a protein) and is asked to provide the totality of concepts with which this object might be associated (i.e., a consistent subgraph). However, the structured-output formulation generally falls under the extreme classification umbrella because the size of the output space is often exceedingly large. This poses problems in measuring similarity between annotations, evaluating accuracy of classification models, and optimization when solving the “argmax problem” <cit.>.

We identify now what we believe is an open problem in computational biology and computer science; that is, efficiently determining the exact number of consistent subgraphs in a given ontology. This problem has a linear-time solution for rooted trees <cit.>, but to our knowledge no such algorithm exists for directed acyclic graphs. This paper therefore proposes a practical solution to this enumeration problem, proves its correctness, analyzes run-time complexity, and introduces various computational speedups. Using this new approach, we analyze four often-used ontologies from the biomedical domain and explore the space of possible annotations. We believe that the algorithms, software, and analysis carried out in this work will lead to better insights into concept annotation spaces and facilitate ontology quality assurance.



§ A MOTIVATING EXAMPLE

A growing number of concept annotation problems are formulated as the manual or computational assignment of a set of mutually related textual descriptors to some objects of interest. One of such problems is the computational prediction of protein function <cit.>, which can be broadly operationalized as follows: 




 
Given: (1) an amino acid sequence with auxiliary data such as structure, expression, interactions, etc. of a protein p with unknown or incomplete function; (2) training data that includes sequences, structures, or systems data corresponding to a (large) set of proteins, some of which have their true biological functions available; (3) a Gene Ontology (GO); i.e., a concept hierarchy used to represent biological functions of proteins in a structured and easy-to-compute-on form.




Objective: provide a set of GO terms that are most likely to be the true (experimental) annotation of p. 
 




The objects of interest here are proteins and the set of textual descriptors of protein function is given by GO – an ontology with a directed acyclic graph structure where each node represents a textual descriptor and each edge represents a particular type of a relational tie between two descriptors <cit.>. 



An example of such an annotation is shown in Figure <ref>, where 8 terms from the molecular function domain have been assigned to this protein. Due to the hierarchical organization of GO, both the set of experimentally determined terms and the set of computationally predicted terms must respect this hierarchy. As shown in this example, the annotation of the term “DNA binding”, implies the annotation of all the other GO terms that conceptually generalize it; e.g., “nucleic acid binding”, “binding”, etc. Typically, the ontology used to represent the annotation space of proteins contains thousands to tens of thousands of terms, whereas the true annotation of a protein consists of tens to at most hundreds of terms. Because the task of a prediction algorithm is to find the most likely annotation, it must devise an efficient procedure to search through the space of all possible annotations.

Most biomedical ontologies have grown over the years to contain a large number of terms. Computationally selecting one such “winning” annotation; i.e., a set of terms, or even providing a short list of most likely annotations, is a significant challenge <cit.>. This prediction problem thus belongs to a so-called extreme classification scenario because the number of possible (discrete) annotations the algorithm must consider is astronomically large. In fact, we noticed that it is not even possible to give an exact number of possible annotations for a protein. Therefore, an answer to such a simple question (“What is the the number of possible GO annotations a protein can be assigned?”) requires the development of a practical counting algorithm. The resulting counts can, in turn, give insight into the nature and the difficulty of the computational function annotation of biological macromolecules.[Reasonable approximations can be provided by calculating the lower and upper bounds, as we have done later in Section <ref>. Neither of those, however, provides a full intellectual satisfaction when an exact count can actually be computed.]

It is important to mention that the annotation of biological macromolecules is one of the most interesting examples of concept annotation, primarily because of its biomedical significance but also because of the sizes of the available ontologies. Similar situations, however, arise beyond computational biology, as in the fields of text mining <cit.> and computer vision <cit.>. 



§ PRELIMINARIES




 §.§ Basic Concepts and Notation

Let 𝒢 = (V, E) be a directed graph, where V is a set of vertices representing concepts and E⊆ V× V is a collection of ordered pairs (u, v) representing directional relationships, u → v, between two concepts. A sequence of vertices u_1, u_2, …, u_k is called a walk if (u_i, u_i+1)∈ E for i = 1, 2, …, k-1. A walk of distinct vertices except for the identical starting and ending vertices is called a cycle. A directed graph that does not contain cycles is referred to as directed acyclic graph (DAG).

Given two vertices u, v ∈ V in a DAG, u is said to be an ancestor of v and v is said to be a descendant of u if there exists a walk from u to v. We denote a set of all ancestors of v as v and a set of all descendants of u as u. We next define v = {v}∪v as the set of extended ancestors of v and u = {u}∪u as the set of extended descendants of u. Finally, if (u, v) ∈ E, the vertex u is said to be a parent of v, whereas v is said to be a child of u. We denote the set of all parents of v as 𝒫(v) and the set of all children of u as 𝒞(u). 



 §.§ Transitivity of Relational Ties

When an object is annotated with ontological concepts, it is often considered that all ancestors of those annotated concepts should be automatically assigned to the object. For example, annotating the function of a protein with “enzyme binding” also implicitly annotates it with “protein binding”, “binding” and, finally, the root term “molecular function”. This type of reasoning requires all involved relationships between concepts to be transitive.

Biomedical ontologies, however, usually contain various types of relationships between concepts, some of which are not transitive. Therefore, we only consider  and  relationships, both of which maintain transitivity and permit reasoning about ancestral concepts. It is also worth noting that we define the direction of edges to be pointing from the general terms to specific so that the depth of a node aligns with the increasing resolution of the descriptors. We show in Section <ref> that the directionality of edges has no impact on the total count. Throughout this work, we consider an ontology 𝒪 = (V, E) to be a DAG, where edges represent transitive relationships.



 §.§ Consistent Subgraphs

Let 𝒪=(V, E) be an ontology and S ⊆ V a set of vertices. A subgraph (S, E^S) is said to be induced from the original graph 𝒪 by S if E^S is the largest subset of pairs (u, v) from E such that both u, v ∈ S. We denote such vertex-induced subgraph as 𝒪S. We also use 𝒪[-S] = (V-S, E^V-S) to denote the subgraph induced by vertices other than S.


A subgraph 𝒪S = (S, E^S) with respect to the original graph 𝒪 = (V, E) is called consistent if ∀ v ∈ S, (u, v) ∈ E  u ∈ S.




§ BASIC ALGORITHMS




 §.§ Problem Specification

Given an ontology 𝒪 = (V, E), our goal is to develop a practical algorithm that enumerates all consistent subgraphs of 𝒪. We allow the graph to have more than a single root (a vertex with no incoming edges) as well as to be disconnected. 

An example of the enumeration problem is shown in Figure <ref>. We generally observe that the number of consistent subgraphs is bounded from below by 2^ℓ, where ℓ is the total number of leaf vertices (those with no outgoing edges), and from above by 2^|V|. The structure of the graph, however, determines the exact count and its proximity to either of the bounds. If the input graph is a chain of |V| vertices (ℓ=1), the total number of consistent subgraphs equals |V|+1. On the other hand, if the original graph is a set of |V| disconnected vertices (ℓ=|V|), there are 2^|V|=2^ℓ consistent subgraphs. This analysis suggests that enumerating consistent subgraphs has a straightforward intractable solution of listing all 2^|V| vertex-induced subgraphs of the ontology and checking for the consistency of each such subgraph.



We use (𝒪) to denote the desired function that takes a directed acyclic graph 𝒪 as input and returns the number of consistent subgraphs in that graph. We use (𝒯) and (ℱ) for the special cases where the input graph is a rooted tree 𝒯 or a forest ℱ, respectively.



 §.§ Counting Consistent Subgraphs in Trees


We first discuss a special case where the input graph is a rooted tree; that is, when each non-root vertex has a single parent. In this case, there exists a linear algorithm in the number of vertices; see Lemma 1 in <cit.>. We provide this solution in Algorithm <ref> with a minor modification resulting from the fact that our algorithm includes an empty tree in the total count. This algorithm naturally extends to collections of rooted trees. One can enumerate subtrees for each tree and take the product as the total count. We refer to this extended algorithm as  (not shown).



Algorithm <ref> recursively traverses a tree in a pre-order manner. For any subtree rooted at vertex v, the number of consistent subtrees that contain v  equals the product of all subcounts from its subtrees rooted at each child. Additionally, we add 1 for the only consistent subtree that does not contain v; i.e., the empty tree. The recursion terminates at the empty tree whose count is one. 





 §.§ Counting Consistent Subgraphs in Directed Acyclic Graphs


Directed acyclic graphs generalize trees in that they allow for multi-parent vertices. 
Such vertices, however, break Algorithm <ref> because the recursive branches are no longer independent. Algorithm <ref> circumvents this problem by recursively decomposing a graph into two strictly smaller subgraphs according to a selected pivot vertex. We will show in the next section that the number of consistent subgraphs in the two smaller graphs add up to be the number for the original graph (Line <ref>, Algorithm <ref>). The algorithm continues recursive enumeration until the graph becomes a forest, in which case it calls . Figure <ref> illustrates the process of graph decomposition with respect to the pivot vertex u. We note that any vertex can serve as pivot and will discuss the selection of pivots and how they impact the run time in Sections <ref> and <ref>.






 §.§ Correctness and Complexity of the Algorithm


We first observe that the size of the problem in the number of vertices is guaranteed to decrease during recursive calls, thus ensuring that the algorithm terminates after a finite number of iterations. Next, we justify the equation corresponding to the Line <ref> in Algorithm <ref>,

    (𝒪) = (𝒪-u) + (𝒪-u).


 
Let (𝒪| u) be the number of consistent subgraphs in 𝒪 that do not contain u. We have
(𝒪| u) = (𝒪-u).



The equal cardinality of the two sets of consistent subgraphs is demonstrated by showing that both sets are contained in each other. For any S ⊆ V-u that induces a consistent subgraph of 𝒪-u, it also induces a unique consistent subgraph in 𝒪. Also, since none of them contains u, we have (𝒪-u) ≤(𝒪| u). Conversely, for any consistent subgraph induced by S such that u ∉ S, we have ∀ v ∈u, v ∉ S by the definition of consistency. Therefore, S also induces a consistent subgraph in 𝒪-u. That is, (𝒪| u) ≤(𝒪-u).



 
Let (𝒪| u) be the number of consistent subgraphs in 𝒪 that contain u. We have
(𝒪| u) = (𝒪-u).



As in Lemma <ref>, for any S ⊆ V-u that induces a consistent subgraph of 𝒪-u, S ∪u also induces a unique consistent subgraph (that contains u) in the original graph. That is, (𝒪-u) ≤(𝒪| u). Also, for any consistent subgraph induced by S and u ∈ S, we have u⊆ S by the definition of consistency. Note that the uniqueness of S implies the uniqueness of S-u. We can see that the subgraph induced by S-u in 𝒪-u is consistent. Given ∀ w ∈ S-u, and (v, w) being an edge in 𝒪-u, we must have v ∈ S-u as well, due to the consistency of 𝒪S with respect to the original graph. That is, (𝒪| u) ≤(𝒪-u).




Given an ontology 𝒪=(V,E) and any u ∈ V, the number of consistent subgraphs in 𝒪 equals the sum of the numbers of consistent subgraphs in 𝒪-u and 𝒪-u.



Equation <ref> holds by combining Lemmas <ref> and <ref>.



To analyze complexity of the algorithm, let n be the number of vertices in the graph and m be the number of multi-parent vertices. Assuming a multi-parent vertex is always selected as pivot, we can express the run time complexity T(n) via the following recurrence

    T(n) ≤ T(n-1) + T(n-3) + f(n),

where f(n) incorporates the time to select the pivot, split the graph and add two large integers. Let us further assume that the larger of the two graphs after decomposition contains n - n/k elements, where 2 ≤ k ≤ n. It is now straightforward to show that T(n) = O(f(n) 2^min(m, s(k))), where s(k)=O(n) if k=O(n) and s(k)=O(log n) if k=O(1).

We can now see that the algorithm is exponential in the worst case; however, it reduces to a polynomial algorithm when m=O(log n) or when k=O(1). Assuming linear time to conduct graph decomposition and a constant time for addition/multiplication, we obtain T(n)=O(n^2). 



§ ADVANCED ALGORITHMS

The run-time of the algorithm heavily depends on the structure of the ontology and the selection of pivots. Here we discuss several practical considerations aimed at accelerating Algorithm <ref>. Once we conclude this discussion, the full method will be presented in  Algorithm <ref> (Section <ref>).



 §.§ Pruning Branching Components

It is easy to observe that when the ontology consists of multiple connected components, these components can be independently and, if needed, simultaneously processed. We take this reasoning a step further to consider a special scenario of nearly disconnected graphs where (i) the two components are connected via a single vertex and (ii) all vertices in one component are descendants of this vertex.


Given a graph 𝒪 = (V,E) and u ∈ V, 𝒪u is called a branching component if  ∀ v ∈ V-u and ∀ w ∈u, (v, w) ∉ E. Vertex u is called a branching vertex.




Figure <ref>a gives an example in which u is a branching vertex, since the removal of u disconnects u (i.e., the branching component, 𝒪^br) from the rest of the graph. We refer to the remaining part of the graph as the stem component, 𝒪^st. More generally, Figure <ref>b shows a graph with a component-wise tree structure, where branching vertices serve as hinges of branching component to their corresponding stems. We will use (𝒪^st, 𝒪^br, u) to denote the desired structure.

Given (𝒪^st, 𝒪^br, u), we demonstrate that (𝒪) can be decoupled into two sequential subproblems: (i) (𝒪^br) and (ii) (𝒪^st). We use φ(u) for the subtotal of consistent subgraphs in the branching component 𝒪^br. We also notice that the entire branching component can be pruned once φ(u) is computed, making u a leaf vertex in 𝒪^st. Therefore, we modify the algorithm so as to allow a subtotal count φ(u) for every vertex as if a branching component has been pruned from u. Notice that φ(u) = 1 for all intermediate vertices and original leaves.

With the introduction of φ(u), the recursive equation in Algorithm <ref> becomes

    (𝒯_r) = φ(r) + ∏_u∈𝒞(r)(𝒯_u).



Similarly, Equation <ref>; i.e., Line <ref> in Algorithm <ref>, must be modified to 

    (𝒪^st) = (𝒪^st-u) + φ(u)·(𝒪^st-u),

where φ(u) accounts for the fact that for any consistent subgraph S_i in the pruned 𝒪^br and any consistent subgraph S_j in 𝒪^st-u, 𝒪S_i∪ S_j ∪u is a distinct consistent subgraph in 𝒪. The approach naturally extends to multiple (hierachical) branching components such that we compute the subtotal of consistent subgraphs within each component and agglomerate them in a reversed topological order.

The pruning operation is preferred before each instance of decomposition for two main reasons: (i) it divides the problem into smaller non-overlapping subproblems, while a direct decomposition usually results in substantial overlapping subproblems; (ii) although a full parallelization over components is restricted since stem components have to be computed only after all of their branching components are finished, the unordered components can be computed simultaneously. For example, as in Figure <ref>b, 𝒪^br_1, 𝒪^br_3, 𝒪^br_4 and 𝒪^br_5 can be computed in parallel.



 §.§ Reverse Graphs

Let 𝒪=(V,E^R) be the reverse graph of 𝒪, where E^R = { (u, v) | (v, u) ∈ E }. We show that the number of consistent subgraphs in 𝒪 equals that in 𝒪.

If 𝒪S is a consistent subgraph of 𝒪, 𝒪-S is a consistent subgraph of 𝒪.



We prove this Lemma by contradiction. For ∀ u ∈ V-S and ∀ v ∈u,[We use u and u for ancestors and descendants of u in 𝒪; u = u and u = u.] if v ∉ V-S, then u ∈v⊆ S due to the consistency of 𝒪S. This contradicts u ∈ V-S. Therefore, the assumption v ∉ V-S is false and we have ∀ u ∈ V-S, ∀ v ∈u⊆ V-S. That is, 𝒪-S is consistent.


This Lemma demonstrates that all complementary white vertices in Figure <ref> form consistent subgraphs in the reverse graph.


Given an ontology 𝒪, (𝒪) = (𝒪).



Given Lemma <ref>, we see that the mapping f(𝒪S) = 𝒪-S is a bijection between the two sets of consistent subgraphs. Therefore, the two sets are of equal cardinality.


Theorem <ref> permits graph reversal at any point during the algorithm depending on which of the graphs is more likely to terminate first. For example, we can always choose the one with fewer multi-parent vertices so as to greedily reduce the upper bound of recursive calls. It is worth noting that all the leaves become roots in the reverse graph. Therefore, in the final algorithm that incorporates both pruning and reversing modules, we generalize the algorithm to allow for φ > 1 on roots (branching vertices in the reverse sense) in order to ensure compatibility.

Having φ(r)>1 on a root indicates that all the ancestors of r have been pruned out. For trees (after pruning), we have 𝒪-r = 𝒪r. With Lemma <ref> and Theorem <ref>, we have

    (𝒪| r) = (𝒪-r) = (𝒪r) = (𝒪𝒟^R(r)) = φ(r).

On the other hand, for any consistent subgraph S containing r, S-r induces a consistent subgraph in 𝒪r and vice versa; thus, 

    (𝒪| r) = ∏_u∈𝒞(r)(𝒯_u).

Hence, these two subtotals sum to be the total count and Equation <ref> remains unchanged. However, if a root r with φ(r)>1 is selected to be the pivot, we have the following equation according to Theorem <ref> and Equation <ref>,

    (𝒪) = (𝒪)    = (𝒪-𝒟^R+(r)) + φ(r)·(𝒪-𝒜^R+(r)) 
       =
        (𝒪-r) + φ(r)·(𝒪-r),

whereas Equation <ref> remains unchanged for non-root vertices.




 §.§ Pivot Selection



As alluded to before, the selection of vertices used for partitioning has the potential to significantly change the computation time. It is therefore reasonable to devise a strategy for pivot selection. Besides a random selection of multi-parent vertices (mpv's), which aims at directly converting DAGs into trees one step at a time, we also consider three other pivot heuristics. The first strategy is to pick a vertex with the maximum degree, with random selection in case of ties, because decomposing the graph according to such vertices may increase the chance of having either disconnected components or branching components. The second strategy selects the pivot so as to minimize e-n+r over the two subproblems, where e, n, and r are the number of edges, vertices, and roots in the two components. We refer to this quantity as “bound” since it is an upper bound of the number of mpv's in the graph (see Supplementary Materials for the proof). Note that it is closely related to the cyclomatic number of the graph. Finally, the third strategy simulates a unit network flow for all vertices running in the direction from leaves to the roots and selects the “bottleneck” vertex; i.e., the one that maximizes the ratio of the flow in the vertex and the number of its descendants (see Supplementary Materials for this pivot selection algorithm). These strategies will be empirically compared in Section <ref>.




 §.§ Hashing

It can occur during the recursive procedure that certain subgraphs require repeated enumeration. In Figure <ref>, for example, the subgraph h-i-j is present in both subproblems shown in Figures <ref>b-c. Computing the count for this subgraph would emerge in the Figure <ref>b subproblem if the ensuing decomposition were based on vertex d, although it would not emerge if the partitioning were based on vertex j. Interestingly, the subgraph k-l would be counted twice in the Figure <ref>b subproblem; i.e., when both d and d are removed, and it would then appear one more time in the Figure <ref>c subproblem.

To avoid repeated enumeration, whenever a solution to a subproblem is obtained, the count for this subproblem is stored. Then, during the recursive calls, we first check if the result is already available before further calculation. To hash a result, we use the sorted IDs of all vertices in the subgraph as a key. Obviously, this key is unique because it corresponds to a vertex-induced subgraph of 𝒪. For the pruned subgraph, we store the key of the subgraph along with the branching vertex. Whenever the ID of the branching vertex is used to generate a key, the stored key of the corresponding subgraph is appended to the vertex's ID with parentheses around it.






§ EXPERIMENTS AND RESULTS


We empirically evaluate the enumeration procedure from Algorithm <ref> and various practical speedups using randomly generated graphs. We then apply this algorithm to four biomedical ontologies to gain insight into the sizes of their concept annotation spaces.



 §.§ Run-time Evaluation


We generated two sets of graphs to investigate the efficacy of our algorithm. Each set contained 1000 graphs with either 25  or 100 vertices. To construct each graph the vertices were added sequentially, with the proposed in-degree in-deg(v) of the k-th vertex v generated according to a Poisson distribution with parameter λ. This vertex then became a child of min(in-deg(v), k-1) previously generated vertices that were themselves selected uniformly randomly. The parameter λ was selected according to the Γ(2.0, 1.0) prior for each new graph and kept constant until the graph was completed. 

With these two sets of simulated graphs, we ran our algorithm with different modules and pivot selection strategies. In particular, we evaluate pivot selection based on (i) random selection of vertices, (ii) random selection of multi-parent vertices, (iii) the degree criterion, (iv) the bound criterion, and (v) the bottleneck criterion. For each pivoting strategy, we subsequently add the pruning component, then hashing, and finally the graph reversal. The criterion for graph reversal was the number of multi-parent vertices; i.e., a graph will be reversed at any point during the recursive process if the reversed graph contains fewer multi-parent vertices.

We report the average wall-time and average number of recursive calls over the two sets of 1000 graphs (|V|=25 in Table <ref>; |V|=100 in Table <ref>). For the smaller graphs, we also ran a brute-force algorithm that was further convenient to empirically evaluate the correctness of our algorithm. We see that simpler schemes perform better on small graphs where the number of recursive calls per graph has not exceeded a few hundreds. On the other hand, the advanced techniques show tangible benefits on the larger graphs reducing the number of recursive calls and total computation time by orders of magnitude. It is possible to envision other variations that could result in further speedups; e.g., selecting multi-parent pivots with the highest degree. These refinements, however, were beyond the scope of this paper.







 §.§ Consistent Subgraphs in Biomedical Ontologies



We use 02/2017 versions of Gene Ontology (GO) and Human Phenotype Ontology (HPO) as the target ontologies and compute the number of consistent subgraphs in each of them. The algorithm is applied to each of the three domains of GO <cit.>: (i) molecular function ontology (MFO; 10,789 terms) (ii) biological process ontology (BPO; 29,575 terms) and (iii) cellular component ontology (CCO; 4,085 terms). Together with HPO (12,167 terms), these four ontologies are widely used in annotating functional terms of gene products <cit.>. We further define the annotation level for each term in the ontology to be the length of the longest path to the root. Starting from the root term, we add more specific terms level-by-level so as to understand how the potential annotation space grows with increased granularity of functional concepts.

In addition to level-wise full ontologies, we also investigate the “used” ontologies in which each term was retained only if at least one protein in the UniProt-GOA <cit.> and HPO <cit.> databases has been confidently assigned that term (confident annotations include all experimental evidence codes as well as “traceable author statement” and “inferred by curators”). Protein function annotations were extracted from the 02/13/2017 release of the UniProt-GOA database, which contains 64,362 proteins with confident MFO annotations, 84,413 proteins with BPO annotations and 79,630 proteins with CCO annotations. HPO annotations were extracted from the 02/24/2017 release of the HPO database where 6,411 genes with confident annotations were extracted.

Figure <ref> shows the completed counts for both full and used level-wise ontologies. For each ontology, we additionally compute the lower bound (generally the larger of 2^ℓ and 2^r, where r is the number of roots) and estimate the upper bound (we convert a graph into a forest by keeping only one randomly selected incoming edge for each multi-parent vertex and then call ).



Although we were not surprised by the astronomical sizes of concept annotation spaces, it was interesting to quantify them whenever feasible as well as to observe an increasing difference between lower and upper bounds (in the 1000s of orders of magnitude) with the level of the ontology. We also find it interesting that a large number of ontological terms have never been used (see Supplementary Materials). These outcomes raise questions regarding the predictability of ontological annotations as most modern algorithms are asked to provide accurate deep annotations to be deemed useful. However, annotation spaces become exceedingly large almost instantaneously with the depth of the ontology, which presents an immense computational and statistical challenge for any prediction algorithm. We therefore believe that the balance between ontology size/complexity and term granularity might become an important topic for future discussions.



 §.§ Entropy of Concept Annotation Spaces


The ability to enumerate subgraphs in relatively large ontologies presents an opportunity to contrast the space of actual ontological annotations in biological databases with the space of possible ontological annotations. To investigate this, we first computed the entropy of actual annotations at different levels in the ontology,


    H(𝒪_lvl) = - ∑_i P(𝒪_lvlS_i) log_2 P(𝒪_lvlS_i),


where 𝒪_lvl is the truncated ontology as in Section <ref>, 𝒪_lvlS_i corresponds to a distinct consistent subgraph annotation observed at that level and P(𝒪_lvlS_i) is the probability that a protein is assigned annotation 𝒪_lvlS_i. We first enumerated all observed subgraphs from the UniProt-GOA or HPO database truncated to a particular level, calculated their relative frequencies, and then plugged these relative frequencies into the entropy formula above. On the other hand, the maximum entropy was computed as log _2(𝒪_lvl) by assuming equal probability for every possible consistent subgraphs.






Figure <ref> shows the ratio between the two quantities for levels greater than 1, suggesting that the world of protein functions, despite great diversity, has low entropy relative to the possible maximum. Although the currently observed functional annotations are incomplete, noisy and biased <cit.>, this suggests considerable departure from the uniform distribution and implies that an extensive number of possible consistent subgraphs have not been used.



§ RELATED WORK

There exists a body of literature in enumerative combinatorics related to our work. One of the most relevant problems is the enumeration of directed acyclic graphs with n distinct (labeled) nodes <cit.>. The resulting count reflects the size of the structure space of Bayesian networks with n random variables and, surprisingly, also corresponds to the number of matrices in { 0,1} ^n× n with all eigenvalues real and positive <cit.>. The number of labeled directed acyclic graphs with n nodes does not have a closed-form solution and is instead available as the A003024 sequence in the On-Line Encyclopedia of Integer Sequences (OEIS); <https://oeis.org/A003024>. The construction was originally proposed by Robinson <cit.> and was further investigated by others <cit.>.

The results on rooted labeled trees include both the enumeration of possible number of trees and also the enumeration of subtrees for a given tree. There are n^n-1 labeled rooted trees with n nodes <cit.> that provide the integer sequence A000169 in OEIS; <https://oeis.org/A000169>. The expansion to forests gives (n+1)^n-1 using Cayley's formula <cit.>, as a single root can be added to connect a forest of unrooted labeled trees into a rooted labeled tree. The recurrence for the number of subtrees of a given tree was proposed by Ruskey <cit.>; see Algorithm <ref>. The generalization to weighted subtrees was given by Yan and Yeh <cit.>. Both algorithms are linear in n assuming constant time addition and multiplication.

Our work also relates to the research in ontology quality assurance. These efforts typically include the analysis of irregularities and redundancy in concept descriptors and graph structure <cit.>. Our work, primarily the software we developed, contributes to this area by facilitating the analysis of the annotation space.



§ CONCLUSIONS

This work presents a practical algorithm for enumerating consistent subgraphs of directed acyclic graphs. Although this study largely addresses an intellectual challenge of efficient subgraph enumeration, it might have practical utility for the studies of annotation spaces in the biomedical sciences and beyond; e.g., in text mining <cit.> and computer vision <cit.>. As ontologies are easy to grow and hard to manually interrogate, we provide a practical tool that can give insights into the complexity of concept annotation tasks <cit.>. As such, it may serve as a guide to ontology developers.



§ ACKNOWLEDGEMENTS

This work has been supported by the National Science Foundation grant DBI-1458477 and the Indiana University Precision Health Initiative. 

abbrv


