Learning directed acyclic graphs from large-scale genomics data.

In this paper, we consider the problem of learning the genetic interaction map, i.e., the topology of a directed acyclic graph (DAG) of genetic interactions from noisy double-knockout (DK) data. Based on a set of well-established biological interaction models, we detect and classify the interactions between genes. We propose a novel linear integer optimization program called the Genetic-Interactions-Detector (GENIE) to identify the complex biological dependencies among genes and to compute the DAG topology that matches the DK measurements best. Furthermore, we extend the GENIE program by incorporating genetic interaction profile (GI-profile) data to further enhance the detection performance. In addition, we propose a sequential scalability technique for large sets of genes under study, in order to provide statistically significant results for real measurement data. Finally, we show via numeric simulations that the GENIE program and the GI-profile data extended GENIE (GI-GENIE) program clearly outperform the conventional techniques and present real data results for our proposed sequential scalability technique.

Introduction

Genetic interaction analysis aims at uncovering the interactions among a set of genes with respect to a specified cell function of a biological system, e.g., the fitness of a specific bacteria colony. The interactions among the genes under study can be characterized by a directed acyclic graph (DAG) [1] where the hierarchical relationship among two genes of a DAG describes their hierarchical interaction type [2]. However, DAGs cannot be observed directly but only the specified cell function under study which yields observable phenotypes. The term phenotype generally describes the specific manifestation of a biological attribute of an organism which can be observed, e.g., for bacteria, a common biological attribute is the growth measured in colony size, where a specific size of the bacteria colony is a phenotype of this biological attribute.

The role of the studied genes in the cell machinery and the hierarchical interaction types of the genes, as well as the DAG, which describes the latter ones, can only be learned by means of knockout experiments where a gene or a set of genes is functionally switched off and the phenotype is observed. Traditionally, only single-knockout (SK) experiments have been conducted but those mainly provide evidence on the importance of a single gene for the investigated cell process and do not convey much information about the interaction among the genes under study.

Recently, with the technological advances in microarrays and the development of the synthetic genetic array technologies [3], new approaches have been taken that are based on large-scale knockout experiments of pairs of genes. Such double-knockout (DK) experiments are much more powerful for exploring genetic interactions since a DK phenotype of an arbitrary pair of genes generally differs considerably from the superposition of the corresponding SK phenotypes of this pair of genes. According to [2], the gene pairs can be classified into one out of five hierarchical relationship classes based on their SK and DK phenotypes. Further, based on the hierarchical relationship classes, the DAG underlying the observed SK and DK phenotypes can be inferred which directly reflects the genetic interactions among the genes.

In order to detect the DAG underlying the SK and DK phenotypes, a variety of statistical methods based on scoring the measurements or on thresholding the genetic interaction (GI)-profile data, which is commonly based on Pearson correlation of the SK and DK phenotypes [4-9], respectively, have been developed. However, methods as presented in [4-9] have three considerable disadvantages: (D1) they show poor performance in detecting the DAG underlying the observed SK and DK phenotypes; (D2) they have no ability to combine different types of side information, e.g., GI-profile data with SK and DK phenotypes, to enhance the detection quality; and (D3) they cannot make use of prior knowledge in order to enhance the DAG detection quality. Especially, the ability to overcome the disadvantage in (D2) will become more important in the future, since there is a constantly increasing amount of different data types, e.g., SK and DK phenotypes, Pearson correlation-based GI-profile data, and other types of GI-profile data, freely available. Furthermore, the ability to overcome the deficit in (D3), i.e., to incorporate a priori knowledge about the existing results in genomics research into the DAG detection procedure, is also of great significance, since existing functional relationships among genes are increasingly better understood based on a variety of studies that constantly extend the knowledge on the cell machinery and molecular biology. Although exhibiting the abovementioned disadvantages (D1) to (D3), methods as those presented in [4-9] are the most commonly used methods to detect the DAG underlying the measured SK and DK data. Therefore, we propose the Genetic-Interactions-Detector (GENIE) program, that is an approach based on the biological system model of [2] with which it is possible to overcome the abovementioned shortcomings of the most popular methods as those reported in [4-9]. Since the hierarchical relationship classes are mutually dependent, classifying each pair of genes to a specific hierarchical relationship class corresponds to a multi-hypothesis test. Thus, we formulate this multi-hypothesis test as a linear integer optimization program [10-15] in order to find the set of hierarchical relationship classes, best matching the observed SK and DK phenotypes. Based on the detected set of hierarchical relationship classes, the set of edges of the DAG which reflects the interactions among the genes can be computed. Furthermore, we propose the GI-GENIE program where we advance the proposed GENIE program by incorporating GI-profile data, e.g., GI-profile data based on Pearson correlation of the observed SK and DK phenotypes, into the DAG detection procedure. Due to incomplete knowledge about the true dependencies among the very most sets of genes, i.e., the true DAG of a set of genes with respect to a specific cell function is unknown or only partially known for almost all sets of genes irrespectively of the cell function under study, there is a strong interest in the genomics research community in statistically reliable statements about the topology of the DAGs underlying large sets of genes, i.e., for the empirical probability of a pair of genes to interact with each other. Towards this aim, we propose a sequential technique based on the GENIE/GI-GENIE algorithms that yields statistically significant statements about the interactions among genes from a large set of genes under study.

This paper is organized as follows. We first summarize the biological system model of [2] in Section 2, and then, we present in Section 3 the GENIE program for detecting the set of hierarchical relationship classes, which represents a valid DAG and matches the DK measurements best. In Section 4, we extend the GENIE program with GI-profile data (GI-GENIE). In Section 5, we present our scalability approach in order to obtain statistically significant results for large sets of genes. Following Section 5, we present results for simulated data which demonstrate the performance of the GENIE and the GI-GENIE methods in Section 6. Furthermore, in Section 6, we display real data results for the scalability approach described in Section 5. Finally, we summarize in Section 7 the key parts of this paper and give a brief outlook on future work.

System model

In this section, we provide a mathematical description of a DAG as well as its biological implications. Furthermore, we introduce the common biological terms and provide a compact description of the genetic interaction model of [2] including simple explanations on how to read and interpret a DAG.

Graph properties of a DAG

According to [16], a graph is well defined by a set of nodes and a set of edges where {a,a,} for denotes a directed edge from a to a and cardinality denotes the number of elements of set . The operators V(·) and E(·) applied to graph yield the set of nodes and the set of edges respectively. We mostly address sets and by and , respectively, for the sake of notational convenience, i.e., . Assume that there is a path P from node to node in graph , i.e., a directed connection from node to node . Then, path P is described by the concatenation of nodes being passed through on the way from node to node , i.e., P=a …a and V(P)={a ,…,a } denotes the set of nodes of path P [16].

The functional dependencies among a set of genes , with elements, for a given cell process and specie can be characterized by a genetic interaction map (GI map,[17-20]) which is essentially a DAG with a common root node, i.e., the reporter level R, [21]. In particular, an arbitrary DAG can be described as a graph with the set of nodes and the set of directed edges . As the genetic interactions can only be observed through the reporter, all edges are always orientated in such a way that each path parting from any arbitrary gene always terminates in the root node R and any gene appears on the path at most once, i.e., there exist no cycles in the graph. Hence, the DAG is always connected via its root node R. For the sake of notational convenience, in most cases, we write gene i when addressing gene g , [21]. The reporter node R is an artificial node, i.e., not a gene, in the concept of a DAG and represents the measured phenotype of the specific cell process under study.

To provide a better understanding of the information encoded in a DAG, we state a simple example, which is similar to the one in [21], based on DAG displayed in Fig. 1. In , there exists an direct edge from gene i 0 to gene j 0, i.e., , which indicates that the activity of gene i 0 controls the activity of gene j 0. Hence, gene i 0 only affects the phenotype via gene j 0 and not directly. We emphasize that in this model, the existence of edge {i 0,j 0} in the DAG only describes the hierarchical functional dependency between genes i 0 and j 0 and not the quantitative effect of gene i 0 on gene j 0.

Fig. 1

DAG of 13 genes and root node R

Biological interaction model

Let us denote as the phenotype for a single gene functionally disabled. In the same way, we define the phenotype for the DK of genes as . Let the datasets and contain all DK phenotypes involving genes . The GI-profile data ρ(i,j) for genes can be computed as the Pearson correlation between the samples of the datasets and , respectively. We remark that the GI-profile data ρ(i,j) does not have to be separately computed as the Pearson correlation of and , respectively. It is commonly extracted from a database where a priori knowledge about the set of genes under study, i.e., , is stored. Since the gene pairs i,j and j,i are identical, it is sufficient to consider only gene pairs . Throughout this paper, we mostly omit the specification that j is greater than i for notational convenience. In genomics research, it is a common assumption that if there is an edge between two genes i,j in DAG , i.e., there is an interaction between genes i,j in DAG , then the GI-profile ρ(i,j) is very likely to be high. Furthermore, according to [2], we assume that each pair of genes i,j belongs to exactly one out of five hierarchical relationship classes that are characterized in Fig. 2. The hierarchical relationship classes are defined according to the model μ (i,j) in which the single-knockout phenotypes R(i) and R(j) are related with the DK phenotype R(i,j). If the gene pair i,j belongs to the hierarchical relationship class k, then the observed DK phenotype R(i,j) is described by the model μ (i,j) provided in Fig. 2. We remark that the five hierarchical dependency graphs in Fig. 2 do not reflect the absolute adjacency relations, but the hierarchical relations between genes i,j in DAG . Hence, given that two genes i,j of DAG are in class k, we cannot conclude that genes i,j are directly arranged in DAG as displayed by the depiction of class k in Fig. 2. This follows from the fact that the description of the hierarchical relationship classes provided in Fig. 2 only contains relative topology information about two genes i,j in DAG . In the following and in addition to [2], we provide, for clarity of presentation, a formal description of the hierarchical relationship classes depicted in Fig. 2 using a graph theoretical representation. Assume that there are I paths P,forτ∈{1,…,I}, from gene i to the reporter node R in DAG and the set containing all such paths is defined as . Furthermore, set contains all J paths from gene j to the reporter node R in DAG . Given gene pair i,j in DAG , then pair i,j belongs to the hierarchical relationship class if and only if condition C as defined below is satisfied:

Fig. 2

Possible hierarchical relationship classes between two arbitrary genes i,j of DAG according to [2]

As stated in condition C1 in (1a), two genes i,j in DAG belong to the hierarchical relationship class k=1, if all paths from gene i to the reporter node R pass through gene j. Hence, gene j is always an element of the set of nodes of each path from gene i to the reporter node R, i.e., j∈V(P) for all paths P from gene i to the reporter node R. With the same line of argument as used in (1a), two genes i,j in DAG belong to the hierarchical relationship class k=2 if condition C2 in (1b) is satisfied. Two genes i,j in DAG belong to the hierarchical relationship class k=3 and are considered to be independent from each other if condition C3 in (1c) is satisfied which states that there is no path P from gene i to the reporter node R that passes through gene j as well as there is no path from gene j to the reporter node R that passes through gene i. As stated in (1d), two genes i,j in DAG belong to the hierarchical relationship class k=4 if there is at least one path P from gene i to the reporter node R which does not pass through gene j as well as for all paths , there is always a path that is a super-path of the respective . With the same line of argument as used in (1d), two genes i,j in DAG belong to the hierarchical relationship class k=5 if condition C5 in (1e) is satisfied.

Class coupling—example

To illustrate this, let us consider the example DAG of Fig. 1. All paths from gene i 0 to node R pass through gene j 0, i.e., they are in a linear pathway with gene i 0 upwards of gene j 0. Thus, the pair of genes i 0,j 0 belongs to class k=1. Note that with the same line of argument, we conclude that also genes i 0 and l 0 belong to relationship class k=1. Since all paths from gene i 0 to the reporter level R do not pass through gene t 0 and all paths from gene t 0 to the reporter level do not pass through gene i 0, genes i 0 and t 0 belong to the hierarchical relationship class k=3 as given in Fig. 2, which states that genes i 0 and t 0 are independent of each other and the DK phenotype amounts to R(i 0,t 0)=μ 3(i 0,t 0). Finally, let us investigate the hierarchical relation between genes t 0 and n 0 in DAG . Obviously, gene t 0 has (at least) one path to node R which does not pass through gene n 0, i.e., genes only having paths to R that do not pass through gene n 0 do not affect the activity of gene n 0. Since there is (at least) one other path from gene t 0 to R passing through gene n 0, we can reason that genes t 0 and n 0 belong to class k=4. Generally, there are strong implications among the hierarchical relationship classes of [2], i.e., if some pairs belong to a specific class, then this has strong implications for all other pairs. Let us consider the case that DAG was not known and only the hierarchical relationship classes for genes i 0 and j 0, i.e., genes i 0 and j 0 belong to class k=1, as well as the hierarchical relationship class for genes i 0 and g 0, i.e., genes i 0 and g 0 belong to class k=1, were available. By definition of the hierarchical dependency graphs in Fig. 2 and the assumptions that genes i 0 and j 0 belong to class k=1 as well as that genes i 0 and g 0 belong to class k=1, we conclude that all paths from gene i 0 to R pass through genes j 0 and g 0. Thus, either all paths from gene g 0 to R pass through gene j 0 or all paths from gene j 0 to R pass through gene g 0. Consequently, genes j 0 and g 0 either belong to the hierarchical relationship class k=1, or k=2.

As we have emphazised by the example above, generally, if the hierarchical relationship class is known for two arbitrary genes i,j as well as for another pair , then this has strong logical implications on the hierarchical relationship classes genes can belong to. Since we can interpret the classification of the pairs of genes i,j, based on their observed SK and DK phenotypes R(i),R(j) and R(i,j), respectively, to exactly one out of the five hierarchical relationship classes as a coupled multi-hypothesis test, we address this problem in Section 3 by a linear integer optimization program. The proposed linear integer optimization program identifies the most consistent set of hierarchical relationship classes, i.e., the set of hierarchical relationship classes that represents a valid DAG and matches best the DK measurements with respect to the logical coupling between the classes. Furthermore, in Section 4, we extend the GENIE program proposed in Section 3 by incorporating GI-profile data in order to jointly detect the most consistent set of hierarchical relationship classes and the corresponding DAG topology.

GENIE algorithm

In this section, we formulate the problem of classifying the gene pairs i,j into the classes of hierarchical relationships based on the observed SK and DK phenotype values as a linear integer optimization program. Furthermore, we translate the logical implications among the hierarchical relationship classes into constraints that ensure that the detected set of hierarchical relationship classes represents a valid graph. That is, the detected set of hierarchical relationship classes represents a graph which is a DAG as defined in Section 2. Finally, we propose a policy to derive an estimate of the true set of edges of DAG based on the detected set of hierarchical relationship classes.

Hierarchical relationship class detection

In order to quantify the mismatch between the measured DK phenotypes R(i,j) and the phenotype model μ (i,j) of class according to Fig. 2, under the hypothesis that the gene pairs i,j belong to class k given their respective SK values, we propose a simple quadratic score [2, 21], as given in Eq. (2)

Let us define the following class-selection variables1

We remark that every DAG can be represented by a set of hierarchical relationship classes which directly corresponds to a set of class-selection variables . The GENIE algorithm of classifying the gene pairs i,j into the set of hierarchical relationship classes that represents a valid DAG and matches the observed SK and DK phenotypes best can be formulated as

where denotes the solution of program OGENIE in (4) and the set of best matching selection variables corresponds to the most consistent pattern of hierarchical relationship classes. Problem OGENIE in (4) is a linear integer program which can be solved efficiently by BB methods [22-29]. The objective of problem OGENIE is to minimize the overall mismatch in classifying each gene pair i,j to one out of five hierarchical relationship classes. The constraints in (4b) reflect the binary nature of the class-selection variables α (i,j), , while (4c) represents a multiple choice constraint to enforce that the gene pairs i,j are only classified to one out of the five hierarchical relationship classes. The set in (4d) is comprised of additional constraints to ensure that the detected set of selection variables always represents a valid graph, i.e., a DAG. In the following, we exemplarily derive topology constraints in set . In order to identify the numerous logical dependencies among the class-selection variables for all , we proceed in the following way. We first fix the assumption that genes i,j belong to class k=1, i.e., α 1(i,j)=1. Further, we assume that genes belong to class k ′, i.e., . Then, we derive the set of classes that genes can belong to under the assumptions made. In the following, we have formulated the logical dependencies among the selection variables for α 1(i,j)=1, i.e., the case that gene i is linearly upstream of gene j, as linear integer inequalities defined in constraints (5a)–(5e) and summarize them as set

where constraints (5a)–(5e) are linear after the continuous relaxation of the selection variables α (i,j). To explain the origin and the functionality of the constraints in , let us further define a sub-genetic interaction map (SMAP) , [21], as given in Fig. 3 according to the following definition where we adopt the graph notation of [16]:

Fig. 3

Example SMAP

Definition 1

Given a non-empty set of edges and a non-empty set of edges , graph , with set of nodes and set of edges , is a SMAP if the following conditions are fulfilled: (i) the graph is acyclic and directed and (i i) there are such that each path P through graph incides via egde e in and leaves graph via edge e out.

DAG , as displayed in Fig. 4, consists of genes i,j and SMAPs and . Obviously, genes i,j belong to class k=1, i.e., α 1(i,j)=1. Furthermore, all genes for which α 1(i,l)=1 must be either located in SMAP or . Thus, it follows from DAG in Fig. 4 that the gene pair j,l is either in hierarchical relationship class k=1 or k=2, i.e., α 1(j,l)=1 or α 2(j,l)=1.

Fig. 4

Schematically reduced DAGs to corresponding to Eqs. (5a)–(5e), respectively

This logical implication is directly reflected by constraint (5a). Given α 1(i,j)=1 and α 1(i,l)=1, the right-hand side (RHS) of (5a) amounts to 1. In this case also, the left-hand side (LHS) of (5a) becomes 1 to fulfill the inequality (5a). Thus, either α 1(j,l)=1 or α 2(j,l)=1. Reversely, assume that α 1(i,j)=1 and α 1(i,l)=1 does not hold, and then, the RHS of (5a) is less than 1, i.e., 0 or −1, while the LHS of (5a) is always greater than 0. Hence, constraint (5a) is fulfilled irrespectively of the choice of α (j,l), i.e., constraint (5a) enforces no logical implications.

Similarly for DAG in Fig. 4, it is obvious that genes i,j belong to the hierarchical relationship class k=1, i.e., α 1(i,j)=1. All genes which are in a linear pathway upstream of gene i, i.e., α 2(i,l)=1, must be located in SMAP . Hence, it directly follows from DAG that also, gene l must be in a linear pathway upstream of gene j, i.e., α 2(j,l)=1. This logical implication is compactly represented in constraint (5b). Under the assumption that α 1(i,j)=1 and α 2(i,l)=1, the RHS of (5b) amounts to 1 enforcing α 2(j,l)=1, so that the LHS of (5b) equals the RHS and the inequality in (5b) is fulfilled. Reversely, assume that α 2(i,l)=0, then the RHS of (5b) is less than 1 and hence the LHS of (5b) is always bigger than or equal to the RHS irrespectively of the choice of α (j,l), i.e., constraint (5a) enforces no logical implications. We can proceed in the same fashion to explain constraints (5c)–(5e) based on DAGs to as given in Fig. 4, respectively. Note that the DAGs to are sufficient illustrations of Eqs. (5a)– (5e) in order to derive all logical implications for the case that genes i,j are in class 1, i.e., α 1(i,j)=1. In the case of DAG , for instance, there can be other DAGs than DAG indeed where genes i,j are in class 1 and genes i,l are in class 1, i.e. α 1(i,j)=1 and α 1(i,l)=1. However, DAG contains all the necessary information in order to derive the logical implications for gene pair j,l, given that α 1(i,j)=1 and α 1(i,l)=1. The same holds for DAGs to . Furthermore, with the same line of argument, we can derive the sets for which reflect the logical implications among the selection variables under the assumptions that α (i,j)=1 for . However, due to space limitations, we omit the derivation of the full set of logical implications at this point and refer the interested reader to [30] where we will provide the full set of topology constraints as well as further supplementary material. The full set of topology constraints in (4d) can be computed as

Finally, a considerable advantage of the presented algorithm is its ability to incorporate prior knowledge into the classification of the gene pairs i,j to the most consistent hierarchical relationship classes. Assume that it is known from existing experimental results that two specific genes do not interact with each other. Then, we can easily incorporate this prior knowledge into program OGENIE in (4) by adding Eq. (7) as defined below

as a topology constraint to program OGENIE. This property is very valuable since it allows the GENIE algorithm to take advantage of existing results in genetic interaction research to improve the reliability of the classification.

Edge computation

Based on the detected set of selection variables which corresponds to the most consistent pattern of hierarchical relationship classes given the observed SK and DK phenotypes, an estimate of the true set of edges of DAG can be computed. It can be theoretically proven that the representation of an arbitrary DAG by its corresponding set of hierarchical relationship classes is not unique. the set of selection variables which directly corresponds to the hierarchical relationship class pattern of DAG represents not only the true DAG but also a set of similar DAGs which have minorly different sets of edges compared to the true DAG . Assume we are only given that for two arbitrary genes i,j of DAG , then we suffer an information loss on the number of paths from gene i to the reporter node R which are independent of gene j. Similarly, given that for two arbitrary genes of DAG , we suffer an information loss on the number of paths from gene j to the reporter node R which are independent of gene i. Hence, this information loss yields ambiguities in computing the set of edges of DAG based on the . In order to clarify the origin of the ambiguities further, let us turn to a simple example. Given DAG as displayed on the LHS of Fig. 5 and the corresponding set of hierarchical relationship classes represented by the corresponding set of selection variables . Note that due to edge and . Now, assume that we want to compute the topology of DAG , i.e., the set of edges , based on . DAG displayed on the RHS of Fig. 5 shows the estimated topology of DAG based on . It can be shown that the black edges of DAG are necessary such that DAG is represented by . Edges e 1 and e 2 in DAG are optional in a sense that their existence has no effect on set . Edges e 1 and e 2 create two paths from gene g 1 to the reporter node R which are independent of gene g 2 and gene g 3, respectively. However, due to edge e 0, gene g 1 already has a path to the reporter node R which is independent of genes g 2 and g 3. Since do not contain information on the number of paths from gene g 1 to R that are independent of g 2, g 3 and g 4, edges e 1 and e 2 do not affect the pattern of hierarchical relationship classes representing DAG , i.e., , and hence, this yields ambiguities in computing the topology of DAG based on its corresponding set of selection variables . Since it is a common assumption in genomics research that GI maps, i.e., DAGs, are not overly dense but rather sparse, we propose a policy which computes the sparsest DAG topology based on the detected pattern of hierarchical relationship classes. Given the detected pattern of hierarchical relationship classes of a DAG , i.e., , we compute an estimate of the true topology set of DAG according to the policy depicted in Table 1 where we make use of the symmetry properties , , , , and . Note that we redundantly expand the set of detected class-selection variables from all pairs to all pairs in order to obtain a compact formulation of the mutually exclusive conditions E1 to E4 as stated in Table 1.

Fig. 5

Left: Original DAG with corresponding set of hierarchical relationship classes . Right: Reconstruction of DAG based on

Table 1

Proposed sparse edge detection policy

Assume that either condition E1 or condition E2 is fulfilled, then we conclude that there is an edge from gene i to gene j in DAG . Given that either condition E3 or condition E4 is fulfilled, we conclude that there exists an edge from gene j to gene i in DAG . We remark that there cannot be an edge between two genes i,j if they are independent of each other, i.e., .

As described by E1, there is an edge from gene i to gene j in DAG , if gene i is linearly upstream of gene j, i.e., , and there is no gene l in DAG that is linearly downstream of gene i, i.e., , and linearly upstream of gene j, i.e., . According to condition E2, there is an edge from gene i to gene j in DAG , if gene i is upstream of gene j with at least one path from gene i to R which is independent of gene j, and furthermore, there is no gene l in DAG that is either linearly downstream of gene i or downstream of gene i with gene i having at least one path to R that is independent of l and neither gene l is linearly upstream of gene j nor gene l is upstream of gene j with an independent path to R. In order to elucidate the effect of condition E2 onto the edge computation, we briefly turn to DAG in Fig. 5. Condition E2 ensures that the optional edges e 1 and e 2 are not detected but only the necessary edges displayed in black color. We remark that conditions E3 and E4 can be elucidated by the same line of argument as used for conditions E1 and E2, but due to space limitations, we omit a detailed explanation at this point. Finally, we propose a condition from which all reporter node edges, i.e, all edges that connect gene with reporter node R in DAG , can be computed. Based on the detected set of hierarchical relationship classes, i.e., , we follow our policy of computing the necessary edges only. For clarity of presentation and notational compactness, we define set as

containing all genes which are in class k=4 with gene , i.e., . Furthermore, we define set as

which contains all genes l of set that are independent of at least one other gene of set . Based on sets and , we formulate condition E as stated in Table 2. We conclude that there is an edge from gene i to reporter node R in DAG , if condition E is fulfilled. Given that gene i is linearly upstream of at least a single gene l, i.e., , there cannot exist an edge from gene i to reporter node R in DAG , since all paths from gene i to R pass through at least one other gene l. Conversely, if there is no such gene l that , then the LHS of E as given in Table 2 is fulfilled. The RHS of E accounts for our policy of detecting sparse DAGs only and is fulfilled if either , , or and are empty. Note that given , it follows that as well, whereas the opposite is not true. In order to explain the effect of the RHS of condition E in an intuitive manner, let us turn to DAG as displayed in Fig. 6. Assume that we are given the pattern of hierarchical relationship classes that corresponds to DAG , i.e., , and we want to compute all reporter node edges based on , i.e., all edges that directly connect a gene in DAG with the reporter node R. Note that , , , and . It can be shown that gene g 3 fulfills the LHS of E, i.e., there is no gene which is linearly downstream of g 3, and furthermore, and . Hence, condition E is fulfilled and edge e connecting g 3 and R is computed in DAG that is the reconstruction of DAG based on . Furthermore, for set , the edge e in the reconstructed DAG is necessary, since , , and . In contrast to e , edge e is not necessary for to represent , since irrespectively of edge e . Hence, e is not detected, since and .

Table 2

Proposed reporter node edge detection policy

Fig. 6

Example DAG to elucidate the functionality of the RHS of condition E of Table 2

We obtain an estimate of the true set of edges of DAG by setting and evaluating conditions E1 to E4 and condition E as stated in Tables 1 and 2, respectively.

GI-GENIE algorithm

In this section, we present the proposed GI-GENIE algorithm which jointly formulates the gene pair classification and the corresponding DAG topology estimation. Let us define the following edge-selection variables

where β(i,j)=1 denotes that there is an edge between genes i,j in DAG and β(i,j)=0 denotes that there exists no edge between genes i and j. Note that unlike α (i,j)=1 for , β(i,j)=1 does not capture directionality information about the graph topology, i.e., β(i,j)=1 states that there is an edge between genes i,j in DAG without specifying the hierarchy among both genes. The topology of any DAG can be represented by the corresponding set of class-selection variables together with the corresponding set of undirected edges {β(i,j)} for all . The set {β(i,j)} can be viewed as the undirected “skeleton” of the DAG that is represented by its corresponding set of class-selection variables . The GI-GENIE algorithm yields an estimate of the true DAG topology by computing sets and which are estimates of the true set of class-selection variables and edge-selection variables, and {β(i,j)}, respectively. Based on SK, DK, and GI-profile data, the proposed GI-GENIE algorithm is formulated as the following LIP:

where the scalars λ ,λ ,λ , and λ are non-negative weighting constants to balance the impact of the SK and DK measurements and the GI-profile data, respectively, on the estimates. In particular, the parameter λ is used for dual purpose: (i) to scale the domain of the knockout scores s (i,j) to the range [0,…,1] which is comparable to range of the correlation data ρ(i,j) and (ii) to trade-off the impact of the knockout scores s (i,j) on the estimation outcome. The parameters λ and λ are in the interval [0,1] where λ ≥λ . The GI-profile (GIP) term is given by

The quotient of defines the threshold for reward of the GI-profile (GIP) term in Eq. (12), where setting the edge selection variable β(i,j)=1 is rewarded if the corresponding GI-profile measurement ρ(i,j) is above the quotient .

The auxiliary variables are generally necessary to ensure that the information about the topology of DAG , which is encoded in the pattern of selection variables detected by program OGI-GENIE, is not contradicting with the set of edge selection variables detected by program OGI-GENIE. In particular, given that the detected pattern of selection variables enforces that there is an edge between genes i,j in DAG , then the auxiliary variables ensure that the corresponding edge selection variable indicates that there is an edge between genes i,j, i.e., . Furthermore, given that the detected pattern of selection variables enforces that there is no edge between genes i,j in DAG , then the auxiliary variables ensure that the corresponding edge selection variable indicates that there is no edge between genes i,j, i.e., . On the contrary, assume that the detected edge selection variables enforce that there is an edge between genes i,j in DAG , i.e., , then the z (i,j) ensure that the detected pattern of selection variables must fulfill one of the conditions stated in Table 1. Consequently, given that the detected edge selection variables enforce that there is no edge between genes i,j in DAG , i.e., , then the z (i,j) ensure that the detected pattern of selection variables does not fulfill any of the conditions stated in Table 1.

Furthermore, let the auxiliary parameters

describe the detection of the edges of DAG based on GI-profile data only, where q(i,j)=1 denotes that there is an edge between genes i,j and q(i,j)=0 denotes that there is no edge between genes i,j. Since any pattern of hierarchical relationship classes implies a specific set of edges and any set of edges implies a specific pattern of hierarchical relationship classes, there is a strong coupling of constraints, i.e., there are strong logical implications among the selection variables α (i,j) and the selection variables β(i,j); the constraints in Eqs. (11e) to (11g) ensure that the detected hierarchical relationship classes and the corresponding edges, i.e., the α (i,j) and β(i,j), are not mutually contradicting. Given that two genes i,j in DAG are independent, i.e., α 3(i,j)=1, there cannot exist an edge between those genes in DAG , i.e., β(i,j)=0. This logical implication is reflected by (11e). Set in (11f) and the linear integer inequalities in (11g) model conditions E1 to E4 of our proposed edge detection policy as stated in Table 1. Since we do not want to redundantly expand the set of selection variables α (i,j) to all in order to not increase the complexity of program OGI-GENIE, we have to consider three cases when formulating conditions E1 to E4 of Table 1 as linear integer inequalities, i.e., , and . Then, the constraints in set

model the logical implications among the selection variables , and β(i,j) for . Together with (11g), constraints (14a)–(14b) model condition E1 of our detection policy taking into account the GI-profile information ρ(i,j) via selection variables β(i,j). Assume that based on the SK and DK phenotypes, it is most consistent that α 1(i,j)=α 1(i,l)=α 2(j,l)=1 for at least one gene l in DAG which corresponds to condition E1 being violated. Hence, there cannot exist an edge between genes i and j in DAG . In this case, the RHS of (14a) amounts to 1 which enforces the LHS of (14a) to amount to 1 as well, i.e., β(i,j)=0. Note that for α 1(i,j)=α 1(i,l)=α 2(j,l)=1, (14b) makes no restrictions on z (i,j). Furthermore, assume that for genes i,j, based on the SK and DK phenotypes, it is most consistent that α 1(i,j)=1, but α 1(i,l) and α 2(j,l) are not jointly 1 for all other genes , i.e., α 1(i,l)+α 1(j,l)<2, then there is an edge between genes i,j in DAG according to condition E1. In this case, it is obvious that (14a) is always fulfilled, i.e., there are no restrictions on β(i,j) by (14a). Since α 1(i,j)=1 and α 1(i,l)+α 2(j,l)≤1 for all , constraint (14b) can only be fulfilled if . Hence, this enforces β(i,j)=1 due to constraint (11g). In this case, constraint (14b) forces . Hence, given that , constraint (11g) sets β(i,j)=1.

Given that the GI-profile data strongly supports that there is no edge between genes i,j in DAG , i.e., β(i,j)=0, and α 1(i,j)=1 is most consistent based on the SK and DK phenotypes measured, then it follows from (11g) that there must be at least one for which z (i,j)=0. In this case, with β(i,j)=0, α 1(i,j)=1, and z (i,j)=0, the RHS of (14b) amounts to 1, forcing the LHS of (14b) to amount to 1 as well, i.e., α 1(i,l)=1 and α 2(j,l)=1, which is together with the assumption of α 1(i,j)=1 a combination that violates the existence of a direct edge between genes i and j. Furthermore, note that (14a) does not have any implications on the selection variables α 1(i,j),α 1(i,l), and α 2(j,l) for the case that β(i,j)=0 and z (i,j)=0.

Assume that the GI-profile data strongly supports that there is an edge between genes i,j in DAG , i.e., β(i,j)=1, and α 1(i,j)=1 is most consistent based on the SK and DK phenotypes measured, then according to (14a), there cannot be any gene for which α 1(i,l)=1 and α 2(j,l)=1. Hence, Eq. (14b) can only be fulfilled if . Thus, (11g) is fulfilled with equality. We remark that given α 1(i,j)=1, constraints (14c) to (14i) are always fulfilled, i.e., they do not pose any implications among the selection variables α (i,j) and β(i,j). Together with (11g), the two inequalities in (14c) model condition E3 where we can elucidate their functionality in the same fashion as before. Constraints (14d) to (14g) along with (11g) model a minor modification of condition E2 where we detect not only all necessary edges but also optional edges given that their existence is strongly supported by the GI-profile. Given that the existence of an edge between genes i,j in DAG is not strongly supported by the GI-profile, i.e., q(i,j)=0, constraints (14d) to (14e) along with (11g) model condition E2 which only allows necessary edges to be detected and we can elucidate their functionality in the same fashion as in (14a) to (14b). Note that (14f) to (14g) are always fulfilled for q(i,j)=0, i.e., no implications among the selection variables α (i,j) and β(i,j) are posed. Assuming that the existence of an edge between genes i,j in DAG is strongly supported by the GI-profile, i.e., q(i,j)=1, then the constraints in (14d) and (14e) are always fulfilled, i.e., no implications among the selection variables α (i,j) and β(i,j) are posed by (14d) and (14e). However, constraints (14f) and (14g) pose relaxed logical implications among the selection variables α (i,j) and β(i,j) compared to constraints (14d) to (14e). Hence, given that q(i,j)=1 and α 4(i,j)=1, an edge between genes i,j in DAG is detected if it is allowed by the pattern of hierarchical relationship classes. Constraints (14h) to (14i) along with (11g) model a minor modification of condition E4 where we detect not only all necessary edges but also optional edges given that their existence is strongly supported by the GI-profile. Furthermore, the functionality of constraints (14h) to (14i) can be explained with the same line of argument as used to elucidate constraints (14d) to (14g).

Denote and as the sets of topology constraints that model the logical coupling among the selection variables , and β(i,j) for and and , respectively. Then, the full set of coupled constraints of the selection variables α (i,j) and β(i,j) is given by

where we again refer the interested reader to [30] for a detailed description of . We obtain an estimate of the true set of edges of DAG based on the computed set of edge selection variables of program OGI-GENIE where we infer the directionality of the edges according to . Note that all reporter node edges are computed according to our proposed reporter node edge detection policy as given in Table 2. Since the reporter node is an artificial node in the concept of a DAG, there is no GI-profile data and thus, no edge selection variable according to (10).

Sequential scalability technique

Due to the combinatorial nature of problems OGENIE and OGI-GENIE, the GENIE algorithm and GI-GENIE algorithm, respectively, cannot be applied to the data of large sets of genes , since the number of candidate solutions to problems OGENIE and OGI-GENIE, respectively, grows exponentially with the number of genes. In order to nevertheless obtain statistically significant statements about the interactions among genes in a large set of genes , we propose the sequential scalability (SEQSCA) technique which is based on the GENIE algorithm and the GI-GENIE algorithm, respectively.

In particular, we create a sequence of S subsets of the full set of genes , i.e., , where we estimate the topology of each DAG , underlying the data of the subset of genes , by the GENIE or GI-GENIE algorithm, respectively, in order to translate the estimated topology of DAG into the corresponding adjacency matrix for each s∈{1,…,S}. Based on the computed sequence of adjacency matrices , we iteratively compute the reliability matrix ∈[0,1] of the full set of genes in such a way that each entry describes the empirical probability of an edge to exist between genes , i.e., the empirical probability that genes directly interact with each other, where a value of 0 means that there is an interaction between the considered pair of genes with probability 0 and a value of 1 means that the considered pair of genes interacts with probability 1.

In particular, in each iteration s, we consider a subset of size of the full set of genes , where each gene of is selected from without replacement with equal probability. Based on the selected subset , we compute in each iteration s an estimate of the true topology of DAG , underlying the observed data of the genes in subset , by the GENIE or the GI-GENIE algorithm, respectively. Furthermore, the topology estimate of DAG is translated into the corresponding adjacency matrix . The update of the reliability matrix for iteration s is computed according to Eq. (16)

with ( being the N×N reliability matrix at iteration s, , ∪ κ ={1,…,N } and κ <κ for all iDAGs but to provide statistical statements to which empirical probability two genes interact with each other.

In Table 3, we have summarized the SEQSCA technique. Finally, by means of the SEQSCA technique, we are able to provide statistically significant statements about the interactions among the genes from a large set by using the GENIE or GI-GENIE algorithm, respectively, in a sequential fashion.

Table 3

Summary of the proposed SEQSCA-algorithm

Initialization: M (0)=0 N×N; M_s=0=0_N_S×N_S; frequency counter ni,j(0)=0

Repeat:

1: Select subset G_s of size N S from G; draw each gene from G with equal probability without replacement

2: Update: ni,j(s+1)=ni,j(s)+1 for all i,j∈G_s

3: Estimate the DAG topology E_s of set G_s using GENIE, GI-GENIE, respectively; ⇒M s

4: Update reliability matrix M (s) according to Eq. (16)

7: Update iteration number: s←s+1

Until: s=S;

Set M_i,j=M^(S)_i,j/ni,j(S)∀i,j∈G

Simulation results

In this section, we first demonstrate the performance of the GENIE algorithm and the GI-GENIE algorithm with respect to conventional techniques for simulated data and further provide statistically significant statements on the interactions among the genes from a large set of genes based on real data using the proposed SEQSCA technique. For the implementation of the proposed algorithms, we used the popular CVX interface [31] along with the well-known MOSEK solver [32].

Synthetic data results

We have generated the ideal SK phenotypes for all as well as the ideal DK phenotypes for all according to the model of [2] as displayed in Fig. 2, where we have corrupted the ideal SK and DK phenotypes by independently and identically distributed zero-mean Gaussian noise with variance σ 2. Furthermore, the GI-profile data has been generated on the basis of the SK and DK phenotypes. We compare both the GENIE algorithm and the GI-GENIE algorithm with the well-known GI-profile approach [2, 33], where the Pearson correlation between the GI-profiles of genes i and j is computed and an edge in the DAG is detected if the Pearson correlation is above a pre-defined threshold t corr, where the directionality is inferred from the selection variable α (i,j) corresponding to the least mismatch model μ (i,j). Furthermore, we compare our proposed methods with the solution of program OGENIE without considering set as a constraint, which means simply classifying each pair i,j to the least mismatch scoring hierarchical relationship class based on the SK and DK phenotypes R(i) and R(i,j), respectively, without ensuring that the resulting pattern of hierarchical relationship classes represents a valid DAG.

In order to limit the Monte Carlo simulation time, we consider a total of 10 genes amounting to 225 binary variables and 2670 constraints for the GENIE algorithm and 630 binary variables and 9645 constraints for the GI-GENIE algorithm, respectively. For the GENIE method without considering the consistency constraints in , we have 225 binary variables and 270 constraints. Since we infer the edge orientation for the Pearson correlation-based method from the least mismatch scoring model, i.e., from the GENIE method without considering the consistency constraints in , we have 270 binary variables and 270 constraints.

In Fig. 7, we display the false detection percentage of edges P ed in the detected DAG normalized to the true number of edges as defined in Eq. (17) versus the SNR.

Fig. 7

D ed versus SNR; t corr=0.6; 200 Monte Carlo runs; λ =0.05, λ =1, λ =0.8

In Fig. 8, we display the percentage of undetected edges P mis in the detected DAG normalized to the true number of edges as defined in Eq. (18) versus the SNR, i.e.,

Fig. 8

D mis versus SNR; t corr=0.6; 200 Monte Carlo runs; λ =0.05, λ =1, λ =0.8

Note that in multi-hypothesis testing problems, it is common to view the diagnostic plots in Figs. 7 and 8 jointly to assess the quality of the proposed algorithms. In Fig. 7, we observe that in the low SNR regime, the Pearson correlation-based method performs best in terms of false detection percentage of edges P ed; however, it fails to improve performance with increasing SNR, because for correct directionality information of the edges, this approach relies on the hierarchical relationship classes detected by method OGENIE without considering . Especially in the high SNR regime, the proposed GENIE and GI-GENIE methods clearly outperform program OGENIE without the topology rule set and approach and respectively reach the performance of the Pearson correlation method. However, the very good performance of the Pearson correlation method in terms of false detection percentage of edges P ed according to Eq. (17) comes at the cost of a rather poor performance in terms of the percentage of undetected edges P mis according to Eq. (18) as can be seen in Fig. 8. In particular, in terms of the percentage of undetected edges P mis, all of the proposed methods outperform the Pearson correlation method. Note that in the high SNR regime, the GI-GENIE of combining SK, DK, and GI-profile data yields the best of both worlds, i.e., it shows an equivalent performance as the Pearson correlation method in terms of false detection percentage of edges P ed, as well as an improvement of the strong performance of the GENIE method in terms of the percentage of undetected edges P mis.

Real data results

Since discovering genetic interaction maps, i.e., DAGs, for specific organisms is an ongoing field of research and the knowledge on genetic interactions is far away from being complete, there is generally no ground truth to directly compare with, even not for yeast which is one of the best understood organisms in this aspect. Therefore, we base the evaluation of the detection performance of the GENIE and the GI-GENIE methods on the biological knowledge that genetic interactions are generally rare and furthermore on the successful detection of known interactions provided by the well-known yeast database of [34]. We remark that to obtain statistically significant statements about large sets of genes, we have applied the proposed GENIE and GI-GENIE algorithms, respectively, along with the SEQSCA technique presented above. To demonstrate the benefit of using multiple data types instead of only one data type, we compare the reliability matrix results for SEQSCA and GI-GENIE with SEQSCA and GENIE which only utilizes SK/DK data. We have applied the abovementioned algorithms to the dataset reported in [35] to obtain the reliability matrices for the GENIE-based SEQSCA as well as for the GI-GENIE-based SEQSCA, G and GI, respectively. The phenotypes reported in [35] are colony size measurements normalized to the wild-type size for each particular SK and DK, respectively. Typically, the colony size serves as a proxy for the fitness of the organism under study, which is the actual cell function of interest that cannot be observed. Therefore, the phenotypes in [35] are non-negative real numbers within the range [0,C max], where C max denotes the maximum size dictated by the experiment setup. For computational reasons, we only considered the first 200 genes, i.e., , of the query gene list of [35]. Figure 9 shows G obtained by the GENIE-based SEQSCA. In Fig. 10, we have displayed GI obtained by the GI-GENIE-based SEQSCA. For both results, we decomposed into a sequence of S=5e4 subsets of equal size N =10. In Fig. 9, 78% of the gene pairs i,j considered by G of the GENIE-based SEQSCA interact with each other with an empirical probability of less than 20%, i.e., [ G]≤.2. Hence, the GENIE-based SEQSCA yields approximately sparse results. This is a good performance in terms of sparsity, since it is known from biology that genetic interactions are generally very rare. However, we can clearly see by the reliability matrix GI that the proposed GI-GENIE algorithm predicts genetic interactions with a much lower frequency as compared to the GENIE algorithm, which means a very good performance in terms of sparsity. We have computed the acceptance ratio

Fig. 9

Reliability matrix G; S=5e 4 subsets considered; subset size N =10

Fig. 10

Reliability matrix GI; S=5e 4 subsets considered; subset size N =10; λ =1e3, λ =1, λ =0.85

where N r is the number of interactions found with high significance ([ G],[ GI]≥1−ε) and which are deposited in the database of [34] as well. N t is the total number of highly significant interactions. Given the confirmed interactions at [34] for our set of genes under study, we remark that evaluating the number of confirmed interactions, that we have also found with our proposed method, would not be a reasonable performance metric, since [34] combines knowledge and experimental results of numerous sources. In contrast to that, we only had the dataset of [35] which only considers a particular phenotype, i.e., colony growth. As depicted in Table 4, we have computed Γ for both the GI-GENIE-based SEQSCA and the GENIE-based SEQSCA. It is obvious that the GI-GENIE-based SEQCA outperforms the GENIE-based SEQSCA, since the acceptance ratio for the GI-GENIE-based SEQSCA is significantly higher than the one of the GENIE-based SEQSCA.

Table 4

Acceptance ratios; ε=0.05

Method:	Γ (%)
SEQSCA and GENIE	53
SEQSCA and GI-GENIE	74

Conclusions

In this paper, we have considered the problem of learning the interactions between genes in a genetic network. We have proposed the GENIE algorithm and the GI-GENIE algorithm to reconstruct the DAG underlying the observed data. The GENIE method is purely based on SK and DK data whereas the GI-GENIE method combines SK and DK data with GI-profile data in order to compute an estimate of the true DAG topology. In Section 5, we have presented the SEQSCA technique in order to obtain statistically significant statements about the interactions among a large set of genes under study. Furthermore, we have shown by simulations that the GI-GENIE algorithm outperforms the conventional techniques and the GENIE algorithm due to the combination of multiple data types, i.e., SK/DK and GI-profile data. Finally, based on the SEQSCA technique, we have presented real data results for the GENIE and the GI-GENIE algorithm, respectively, where we have confirmed that the GI-GENIE method outperforms the GENIE method.

Endnote

1 In a discrete optimization context, the class-selection variables defined in (3) are denoted as SOS-1 type variables. However, for the sake of readability, we will mostly omit this optimization context-based annotation and refer to the variables defined in (3) as class-selection variables.