Explaining Genetic Knock-Out Effects Using Cost-Based Abduction

1 Introduction

Since the word gene was coined in 1909, it had been a common belief that the higher an organismÕs complexity, the more genes it has. However, genome sequencing has revealed that the entire human genome contains only 23,000 to 40,000 genes, which is close to the number of genes in some types of worm [Baldi and Hatfield, 2003]. It is now known that it is not the number of genes, but gene interaction and regulation that are the sources of organism complexity. As a result, Genetic Regulatory Networks (GRN) have become one of the most interesting and challenging problems in computational biology and is expected to be the center of attention for a few decades to come [Baldi and Hatfield, 2003].

Other mathematical modeling and approaches to GRN include discrete models like Boolean Networks [Kauffman, 1993], continuous models like ordinary, stochastic and qualitative differential equations [Baldi and Hatfield, 2003], probabilistic and graphical models including Bayesian Networks (BN) [Beer and Tavazole, 2004; Friedman, 2004; Huttenhower et al., 2006] and factor graphs [Yeang et al., 2004] , and rule-based models including Induction Logic Programming [Ong et al., 2007], and Abduction Logic Programming [Papatheodorou, 2007; Ray and Kakas, 2006]

Despite the advantage of being able to integrate different sources of biological knowledge in a single, homogeneous knowledge base, the area of rule-based models has received little attention in computational biology in general and in GRN modeling in particular. This is because there is no direct way to integrate an objective or a cost function with the knowledge base to measure solution quality. This issue is addressed to some extent by probabilistic inductive logic programming, which combines probability and logic programming [Raedt et al., 2008]. In addition, there is no clear mathematical correspondence between rule-based models and other machine-learning approaches, such as probabilistic methods and Neural Networks (NN).

In this paper, we propose using Cost-Based Abduction (CBA) to model GRN by explaining the effects of genetic knock-out experiments. Because CBA is a rule-based system, it integrates different data sources efficiently and easily. In addition, it provides an associated cost for each explanation of the data; this cost serves as an objective function that we want to minimize. Interestingly, CBA has a clear, mathematical correspondence to BN, as we show later. Our results show that CBA is a promising tool in the field of computational biology.

This paper is organized as follows: the rest of the introduction gives the necessary background on CBA, its relation to Bayesian inference and NN, and describes the data sources used in the paper. We then describe how to use the data sources to create an un-annotated, or a skeleton graph. Finally, we describe our algorithm for building a CBA system from the skeleton graph. In particular, we show how to create rules and assign costs in a way that naturally enforces the constraints of genetic pathways and explains the effects of genetic knock-outs.

1.1 Cost-Based Abduction

CBA was first introduced by Charniak et al. [Charniak and Shimony, 1990]. Formally, a CBA system is a 4-tuple

( , , , ) H R c G , where H is a set of hypotheses, or propositions, c is a function from H to a nonnegative real c(h) called the assumability cost of h � H, R is a set of rules of the form: 1 2 :( ... ) n k i i i i R p p p p ∃ ∃ τ for all , .., i n i i p p H ∨ , k i p H ∨ and G H ∧ is the goal, or the evidence set [Abdelbar, 1998].

Our objective is finding the Least Cost Proof (LCP) for the goal. The cost of a proof is the sum of the costs of all the hypotheses that were assumed to complete the proof. Any given hypothesis pi � H can be true either by proving it or by assuming it and paying its assumability cost. Hypotheses that can be assumed have finite assumability costs and are called ÒassumablesÓ. Hypotheses that are proved from the assumables are called ÒprovablesÓ. A provable can be thought of as an assumable with an infinite assumability cost that we cannot afford. We can assume, without loss of generality, that all rule consequents are provables. To see this, suppose i p is an assumable with assumability cost ( ) i c p x _ and suppose it appears as the consequent of at least one rule, i r . We can create a new hypothesis i pζ with an assumability cost equals to x , set the assumability cost of i p to ι and add the rule: i i p p ζ τ . In addition, for the purpose of this paper, we call every hypothesis belonging to the goal set G H ∧ a sub-goal hypothesis. All sub-goals in our system are provables.

Finding the optimal solution in CBA has been shown to be NP-Hard [Charniak and Shimony, 1994]. The most significant approaches to finding LCP for CBA can be found at [Den, 1994; Ishizuka and Matsuo, 2002; Ohsawa and Ohsawa, 1997]. One of the most successful algorithms uses Integer Linear Programming (ILP) as developed by Santos [Santos Jr., 1994]. Abdelbar et al. showed how to solve CBA using High-Order Recurrent Neural Networks (HORN) [Abdelbar et al., 2003; Abdelbar et al., 2005].

Finding the LCP in a CBA system is equivalent to finding the maximum a posteriori assignment (MAP) in BN [Charniak and Shimony, 1990; Charniak and Shimony, 1994]. They showed that everything one can do using BN can also be done using CBA, and vice versa. However, despite their equivalence, CBA has many advantages. For instance, it is believed that finding LCP is more efficient than finding MAP, and it may be easier to find heuristics for CBA than for BN [Abdelbar, 1998; Charniak and Shimony, 1994]. Also, using costs instead of probabilities gives another perspective to the problem and is easier for some people to grasp. We have also found that the CBA knowledge representation in terms of rules is more natural for modeling genetic regulation. Finally, Santos has found necessary and sufficient conditions under which a CBA system is polynomially solvable [Santos Jr. and Santos, 1996]. In contrast, conditions for polynomial solvability of finding MAP in BN have not been found, even with applying restrictions on the graphical representation [Shimony, 1994] and even for trying to find an alternative next-best explanation [Abdelbar and Hedetniemi, 1998]. Andrews et al. provided a framework for finding MAP for BN using HORN and using CBA as an intermediate representation [Andrews and Bonner, 2009]; the theoretical foundation of that framework shows the strong equivalency between BN/CBA and HORN search spaces [Andrews and Bonner, 2011].

1.2 Biological Data Sources

In this paper, we illustrate the application of CBA to GRN by using it to model and learn the well-studied pheromone pathway in Yeast. We use three well-known data sets for this purpose. The first data set consists of protein-DNA interactions, also called location data or factor-binding data, from [Lee et al., 2002]. The directionality of these interactions is known to be from a protein to a gene in the DNA. The second data set consists of protein-protein interactions from the well-known DIP database. The directions of the protein-protein interactions are not known a priori and are learned by CBA. Finally, we use the knock-out data from [Hughes et al., 2000]. Knock-out data describes the effects of induced deletion mutation experiments on some genes. We do not use the knock-out data to build the skeleton graph. Instead, the effect of the knock-outs is what is being explained by the CBA system, and the resultant explanations annotate the entire graph, which in turn represents the learned GRN model. Both protein-DNA and protein-protein data indicate whether there is an interaction without specifying its sign. CBA uses the knock-out data to infer the sign of each interaction, as well as the directions of protein-protein interactions.

2 The Method

Our modeling method, described below, can be summarized as follows. We first select the elements in the genome that we are interested in, which we call the elements of interest set. In this paper, we select all genes and proteins that are known to participate in the pheromone pathway. We then create a skeleton graph. Each node in the skeleton graph represents a gene and its protein product, which is unique across all data sets. An edge in the skeleton graph represents a potential interaction between a pair of nodes. Using this graph, we enumerate all paths between each knocked-out gene and genes it is known to affect. Our CBA builder algorithm then uses these enumerated paths to build a CBA system that effectively annotates the entire graph while looking for the LCP of all the knock-out effects. The resulting annotated graph is the learned model for the GRN under study.

2.1 The Un-annotated (Skeleton) Graph

To build the skeleton graph, we first create a node for each element in our elements-of-interest set. Graph edges are then created in two phases. Phase 1 creates (directed) Protein-DNA edges: we search the location dataset for interactions between nodes of the graph; and for each such interaction, we create a directed edge from the protein to the DNA. Phase 2 creates (undirected) Protein-Protein edges: we search the DIP database for protein-protein interactions between nodes of the graph; and for each such interaction, we create an undirected edge between the two nodes. These two phases result in an un-annotated graph. We then use a recursive breadth-first algorithm to enumerate all possible paths between each knocked-out gene and genes it is known to affect. We call these paths enumerated or potentially valid paths.

The CBA system (described below) adds annotations to the skeleton graph. In the annotated graph, all edges have signs and existence labels, and in addition, undirected edges have directions. As in [Yeang, et al., 2004], a valid path ka Μ starts with a knocked-out gene k g and ends with an observed (or affected) gene a g , and must satisfy the following constraints:

1. All edges are in forward direction, from k g to a g .
2. The aggregate sign of the edges in the path must be consistent with the sign of the knock-out effect.
3. The path must end with a protein-DNA edge.
4. The path must be no longer than a predefined upper bound.
5. If an intermediate gene in the path is also knocked out, it must exhibit a knock-out effect on a g .
6. Each edge must have one and only one direction and one and only one sign.

We use the potentially valid paths enumerated from the skeleton graph to build a CBA system whose LCP enforces these constraints to give us the best valid paths and a fully annotated graph.

2.2 Example

Figure 1 shows a skeleton graph that might be generated by the algorithm described above. bove.

Let us assume that we have two knock-out effects to explain, 1 3 ( , , ) g g ∞ 2 and 3 5 ( , , ) g g ∞ 2 . The notation ( , , ) i j g g ∞ 2 means that deletion of gene i g results in down-regulation of gene j g . We say that a knock-out effect 1 3 ( , , ) g g ∞ 2 is explained if there is a valid path that connects both genes. The skeleton graph above suggests three potentially valid paths that might explain these knock-out effects: 1 2 3 g g g τ τ , 1 4 3 g g g τ τ , and 3 4 5 g g g τ τ . Clearly, the direction of the edge between 3 g and 4 g is vital to any explanation of the data. In particular, the direction has to be from 3 g to 4 g . If this direction is reversed, we will not be able to explain 3 5 ( , , ) g g ∞ 2 . The two paths that explain all the data are 1 2 3 g g g τ τ and 3 4 5 g g g τ τ . In addition, to be consistent with the knock-out effects, the edge signs along both paths must be either ( , ) 2 2 or ( , ) 0 0 . (The aggregate sign of the edges in a valid path must be the opposite of the observed knock-out effect, since the deletion itself is negative.)

3 The CBA Builder Algorithm

This section shows how we use the potentially valid paths to build a CBA system that finds valid paths that provide the best explanation of all the knock-out data and that fully annotate the graph.

3.1 Creating the CBA Hypotheses

For each potentially valid path ka Μ for a knock-out pair ( , , ) k a g g ∞ χ , we create the following hypotheses:

• ♥ kj k 0 , kj k 2 and 0 kj k , for every intermediate gene j g in ka Μ , including the observed gene a g . These hypotheses represent aggregate positive, negative and zero knock-out effects, respectively, from the knocked-out gene k g to the intermediate gene j g . kj k 0 and kj k 2 are provables while 0 kj k is an assumable.
• ♥ ka o . This is a sub-goal hypothesis representing the observed knock-out effect to be explained.

For each edge ij E between i g and j g we create the following hypotheses:

• ♥ 1 ij x and 0 ij x , assumables representing the existence of the edge in the GRN (i.e., whether or not it represents a real biological interaction).
• ♥ ij s 0 and ij s 2 , assumables representing the sign of the edge.
• ♥ i j d τ and j i d τ , assumables representing the direction of the edge. (Used for undirected edges only.)
• ♥ ij s , a sub-goal to determine the sign of the edge.
• ♥ ij d , a sub-goal to determine the direction of the edge. (Used for undirected edges only.)
• ♥ ij y , a sub-goal to determine whether the edge exists in the GRN.

3.2 Creating the CBA Rules

Given a potentially valid path ka Μ for a knock-out pair ( , , ) k a g g ∞ χ , let the genes on the path be 0 1 , ,..., m j j j g g g , where 0 j k g g _ and m j a g g _ . For the first edge 1 kj E on the path, we create two rules:

• ♥ 1 1 1 1 1 kj kj k j kj x s d k 0 τ ∃ ∃ τ

2

• ♥ 1 1 1 1 1 kj kj k j kj x s d k 2 τ 0 ∃ ∃ τ

2

For every subsequent edge 1 n n j j E , we create four rules:

• ♥ 1 1 1 1 1 n n n n n n n n j j j j kj j j kj x s k d k 2 2 2 2 0 0 τ 0 ∃ ∃ ∃ τ
• ♥ 1 1 1 1 1 n n n n n n n n j j j j kj j j kj x s k d k 2 2 2 2 0 2 τ ∃ ∃ ∃ τ
• ♥ 1 1 1 1 1 n n n n n n n n j j j j kj j j kj x s k d k 2 2 2 2 2 0 τ ∃ ∃ ∃ τ

2

2

• ♥ 1 1 1 1 1 n n n n n n n n j j j j kj j j kj x s k d k 2 2 2 2 2 2 τ 0 ∃ ∃ ∃ τ

where the rules contain directionality assumables only if the edge is undirected.

The main feature of these rules is that for any two consecutive edges on the path, the consequents of the rules for the first edge are antecedents of the rules for the second edge. This makes it easy to enforce the second path constraint above, since the consequent of each rule represents the aggregate effect of all the edges along the path from the first edge up to the edge for the rule in question.

For each undirected edge, we create two rules, which intuitively mean that the edge has two possible directions:

• ♥ 1 1 n n n n j j j j d d 2 2 τ τ and 1 1 n n n n j j j j d d 2 2 τ τ

For every edge, we create two rules, which intuitively mean that the edge has two possible signs:

For every edge, we create two rules, which intuitively mean

• ♥ 1 1 n n n n j j j j s s 2 2 0 τ and 1 1 n n n n j j j j s s 2 2 2 τ that the edge does or does not exist in the GRN:

For each knock-out pair ( , , ) k a g g ∞ χ , we create the following three rules, which intuitively mean that the knock out can have three possible effects on the observed gene: upregulation, down-regulation or leaving it unchanged.

• ♥ ka ka ka k ck o 0 0 ∃ τ , ka ka ka k ck o 2 2 ∃ τ , 0 ka ka k o τ

Here, ka ck 0 and ka ck 2 are auxiliary assumables that carry the costs of positive and negative knock-out effects.

Finally, all sub-goal hypotheses 1 1 1 { , , , } n n n n n n j j j j j j ka y d s o 2 2 2 are added to the goal set.

3.3 Assigning Costs to Hypotheses

The last step in creating the CBA system is assigning the proper cost to each assumable. For all ij s 0 , ij s 2 , i j d τ and j i d τ , we assign an equal and relatively high cost, to force the system to assume only one sign and only one direction while not being biased towards a particular choice. We calculate the cost of the existence hypotheses 1 0 { , } ij ij x x and the knock-out effect hypotheses { , } ka ka ck ck 0 2 based on observed experimental values reported in the datasets.

Calculating costs for edge existence assumables

The costs of the existence hypotheses 1 0 { , } ij ij x x for an edge depend on the likelihood that the edge represents a real biological interaction. These likelihoods are derived from the protein-DNA and protein-protein data. The likelihood itself is not reported in the data, but the p-value is.

We recover the likelihood using the same methods as [Yeang, et al., 2004]. For example, for protein-DNA data, we assume the p-value comes from testing the null hypothesis 0 H ( 0 1 ij x _ , the interaction does not occur) against the alternative hypothesis 1 H ( 1 1 ij x _ , the interaction does occur). Accordingly, we recover the log likelihood that the interaction has occurred as follows:

where D is the observed data, 0 d and 1 d are the degrees of freedom in 0 H and 1 H , n is the sample size the p-value was computed from, which is 3 in our case, and L is recovered from the p-value using the following approximation:

where p is the reported p-value, and F is the cumulative 2 ∪ distribution with one degree of freedom. The costs of edge existence are then computed as follows:

Calculating costs for knock-out effect assumables

The likelihood of a knock-out interaction can be recovered using the same approximation mentioned above. If the likelihood of the knock-out effect ka k between genes k g and a g is � , then the cost of assuming no knock-out effect is

The cost of assuming a positive knock-out effect ka k 0 is

The cost of assuming a negative knock-out effect ka k 2 is

Here, ka o is the observed effect of the knock-out on the expression level of gene a g , and ⌡ is a very small real number representing the possibility that the effect is due to something other than the knock-out experiment.

4 Results

We used integer linear programming (ILP) to solve the resultant CBA system. In particular, we used the popular public domain LP-solve to solve the CBA system after converting it to an equivalent ILP program.

In addition to toy examples like Figure 1, we applied our method to the well-studied Yeast pheromone pathway. The skeleton graph consisted of 37 protein-DNA interactions and 30 protein-protein interactions. We tried upper bounds of up to 10 nodes per path. In all cases, the paths discovered to explain the knock-out effects satisfied the valid path constraints mentioned above. Table 1 summarizes the result for different upper bounds:

4.1 Comparison

Many of the elements of our work are based on [Yeang, et al., 2004], which uses a probabilistic approach to infer a GRN of Yeast from multiple biological datasets. They use the same datasets as we do and create what they call a physical network model. The main difference is that they use factor graphs instead of CBA, and they compute MAP instead of LCP. We found that the CBA approach produced better results for the pheromone pathway. For instance, their method was only able to find approximate MAP configurations, while CBA was able to find an exact LCP. Moreover, while they used a maximum path length of 5, CBA could easily handle paths of length 10. Finally, because CBA computes the exact LCP, the solution does not have variant parts as in models that compute approximate MAPs.

A complete discussion and comparison of biological results is beyond the scope of the present paper. However, using CBA, we did recover almost all the confirmed signal transduction directionalities in the pheromone pathway, including STE11τSTE7, {STE5,STE7}τFUS3 and {FUS3,KSS1}τSTE12. The 21 direct regulations of STE12 were also detected. Neither FUS3 nor KSS1 violated any path constraints. We recovered a consistent relation between STE12 and MCM1 in both signals and directionalities.

The next section elaborates on interesting features of the size of the CBA system produced by our build algorithm. These features explain why our model can handle path lengths longer than other probabilistic models.

4.2 CBA Size Features

It is easy to verify that apart from the goal rule, our CBA builder will always produce rules that consist of at most 5 variables, including the consequent. This fixed length guarantees that the complexity of the graph will never cause an unexpectedly long constraint when using ILP or an unexpectedly high-order edge when using NN. This short rule length also makes the job of the CBA solver easier. However, the difficulty of a CBA instance depends on more than just rule length. Other factors include solution depth, the ratio of assumables to provables, and the ratio of the number of hypotheses to the number of rules.

In the probabilistic approach of [Yeang, et al., 2004], the size of the factor graph depends on the number of paths in the GRN. In contrast, in the approach described here, the size of the CBA system depends only on the number of edges in the GRN. This is possible because, as a logical formalism, CBA very naturally deals with both AND and OR constraints and with constraints based on paths and subnetworks in a graph. Consequently, the hypotheses and rules for an edge are created only once, and reused for each path that the edge belongs to. As a result, increasing the number of paths or increasing path length does not result in an exponential increase in the size of the CBA system. Figure 2 illustrates the relation between path length and other aspects of CBA complexity for the pheromone pathway.

5 Concluding Remarks and Future Work

We applied CBA to the modeling of genetic regulatory networks (GRN). CBA can easily integrate biological data from multiple sources and explain the effects of gene knockout experiments. The size of the CBA instance does not increase exponentially with problem size. Our results suggest that CBA can be a useful tool for computational biology in general and for GRN modeling in particular. We tested our algorithm on the pheromone pathway in Yeast. In addition to testing our model on other pathways and larger GRNs, future work includes further optimization of the CBA instance. We need, for instance, to make the path enumeration phase more efficient, either by building the CBA rules directly from the edges without any path enumeration or by using a more efficient algorithm than recursive breadth-first search. It will also be interesting to use the node ordering of the CBA instance to create an equivalent BN and compare performance. Finally, because we use an ILP solver, we can study the characteristics of the resulting constraint matrix and check whether it meets the CBA polynomial-solvability conditions determined by Santos.