Learning a Markov Logic network for supervised gene regulatory network inference

Abstract

Background

Gene regulatory network inference remains a challenging problem in systems biology despite the numerous approaches that have been proposed. When substantial knowledge on a gene regulatory network is already available, supervised network inference is appropriate. Such a method builds a binary classifier able to assign a class (Regulation/No regulation) to an ordered pair of genes. Once learnt, the pairwise classifier can be used to predict new regulations. In this work, we explore the framework of Markov Logic Networks (MLN) that combine features of probabilistic graphical models with the expressivity of first-order logic rules.

Results

We propose to learn a Markov Logic network, e.g. a set of weighted rules that conclude on the predicate “regulates”, starting from a known gene regulatory network involved in the switch proliferation/differentiation of keratinocyte cells, a set of experimental transcriptomic data and various descriptions of genes all encoded into first-order logic. As training data are unbalanced, we use asymmetric bagging to learn a set of MLNs. The prediction of a new regulation can then be obtained by averaging predictions of individual MLNs. As a side contribution, we propose three in silico tests to assess the performance of any pairwise classifier in various network inference tasks on real datasets. A first test consists of measuring the average performance on balanced edge prediction problem; a second one deals with the ability of the classifier, once enhanced by asymmetric bagging, to update a given network. Finally our main result concerns a third test that measures the ability of the method to predict regulations with a new set of genes. As expected, MLN, when provided with only numerical discretized gene expression data, does not perform as well as a pairwise SVM in terms of AUPR. However, when a more complete description of gene properties is provided by heterogeneous sources, MLN achieves the same performance as a black-box model such as a pairwise SVM while providing relevant insights on the predictions.

Conclusions

The numerical studies show that MLN achieves very good predictive performance while opening the door to some interpretability of the decisions. Besides the ability to suggest new regulations, such an approach allows to cross-validate experimental data with existing knowledge.

Background

Gene regulatory network inference has received a lot of attention over the last decade due to the abundance of high-throughput data. A gene regulatory network (see for instance [1]) usually refers to a set of genes whose expression varies over time due to the inhibitive or inductive roles of regulators. Deciphering these regulations at work in the cell will provide a thorough understanding of the cell behaviour and will eventually aid in controlling or repairing when needed. Inference of gene regulatory networks as a problem of empirical inference fits the framework of machine learning as described in [2]. Three main families of inference algorithms have been developed so far: (1) unsupervised model-free approaches that use information theory to extract a non-oriented graph of dependence between variables, (2) unsupervised reverse-modeling approaches that model the network behavior as a (dynamical) system [3] and (3) supervised edge prediction approaches that focus on the graph of regulation and only predict the presence/absence of regulations [4-7]. In the first family, relevance networks like ARACNE [8], CLR [9] and TD-ARACNE [10] use a mutual information score between the expression profiles of each pair of genes and given a threshold, decide to predict an interaction or not. The second family is based on model of behavior of the network, either static or dynamic. In case of static models devoted to steady-state data, Gaussian Graphical Models (GGM) [11, 12] allow to build a linear regression model that expresses how one gene can be predicted using the set of remaining genes. Interestingly, GGM build a network using partial correlation coefficients, providing a stronger measure of dependence compared to correlation coefficients used in relevance networks. A powerful approach to regression and network inference based on an ensemble of randomized regression trees [13] has also proven to outperform competitors in inferring gene regulatory networks in recent DREAM competitions. Bayesian networks [14] provide another important approach in static modeling. Learning a Bayesian network involves learning the acyclic oriented graph that describes the parental relations between variables and the conditional probabilities that govern the behavior of the network. While appropriate to gene regulation cascades, Bayesian networks cannot, however, model cycles in the network. Other models incorporating dynamical modeling have therefore been proposed in the literature: dynamical Bayesian networks and differential equations [15-17].

Taking a different angle, supervised edge prediction methods build a decision function that associates a class label to a pair of vertices (genes or proteins) without searching for a model of the network behavior. These methods assume that the network to infer is partially known and that information on the vertices are available. They have been mainly developed for protein-protein interaction network inference, using kernel methods [18-23]. The principle underlying [20, 21] is to build pairwise Support Vector Machines (SVM) with an appropriate definition of kernels between pairs of proteins from a kernel defined between individual proteins. Pairwise kernels can also be combined into a linear combination (usually an average) to deal with multiple sources of information. In [23], another point of view is taken: local models (still SVMs) are attached to each target protein in order to predict whether a candidate protein interacts with the considered target, and these models are then combined. Recently, the work of [22] has shown that the local model is equivalent to a pairwise SVM considering a local definition of a pairwise kernel.

In the case of gene regulatory network inference, the supervised setting of edge prediction has been explored less. It was first introduced by Qian et al. [4] using gene expression as unique descriptor and further developed by Mordelet et al. with the SIRENE method [5]. Similarly to [23], SIRENE estimates a local model for each transcription factor and then combines all local models together. The method requires a list of known transcription factors that serve as targets. Other advances in supervised edge inference concern with the problem of lack of true negative examples and therefore focus on learning from positive only and unlabeled examples. Some methods develop strategies to select reliable negative examples from the unlabeled set and then solve a classical balanced binary classification problem [24, 25]; others adjust the probability of being positive estimated by a classifier trained on positive and unlabeled examples [6, 7, 26].

Choosing between the three kinds of network inference methods, namely model-free, model-driven and supervised approaches, relies on the goal of the study. Model-free approaches give a good first network approximation when only one kind of data is available. Reverse-modeling delivers a model of the network that can be used to predict its behavior but requires a sufficient amount of observations, if possible acquired with different initial conditions or perturbations. Supervised edge prediction is relevant when a sufficiently large set of regulations is known a priori and various sources of gene annotations are available. It will be especially meaningful when the biologist wants to increase the corpus of existing knowledge.

This paper deals with the latter prediction problem. We assume that a directed graph of regulations is known partially for a target set of genes. For instance, it is the result of the biologists’s experience and careful mining of the literature. Besides the graph structure, we also suppose that a set of various descriptors of genes and their products are available for the target set of genes, such as gene expression data, Gene Ontology (GO) annotation, protein-protein interaction and also genes location on chromosomes. Our goal is to build a decision function that predicts if an ordered pair of regulator and regulee candidates belongs to the class Regulation or No Regulation.

In this work, we address four issues raised by supervised edge prediction and implement the whole approach on a new experimental dataset related to the ID2 genetic regulatory network in human keratinocytes. The first issue concerns the available sources of information about genes and proteins. These sources provide multiple views of the data which are by definition heterogeneous and very often highly structured. The second issue is related to network inference interpretability: many of the proposed methods are black boxes, while biologists are interested in how the predictions have been obtained. The third issue, as raised by many authors, deals with imbalanced data: very few positive examples of “regulation” is available compared to the huge number of negative examples of “no regulation”. Finally, the fourth issue we tackle in this paper, concerns the performance assessment of a supervised edge prediction tool. Although the best performance assessment comes when biologists go back to the experimental laboratory to test prediction of new regulations with additional and independent experiments, there is a lot of room for in silico studies to measure the ability of an edge prediction tool to provide evidence for regulations. The first and the second issue call for a common framework of representation for all the views of the data. For that purpose, we use first-order logic to represent both data and background knowledge. In order to benefit from the tools of statistical learning and to avoid some of the weaknesses of pure inductive logic programming raised, for instance, in [27], we choose a Markov Logic network (MLN) [28, 29] as the edge predictor. MLN allows to make predictions using a set of weighted first-order logic rules, thus providing interesting insights on decisions. The third issue is systematically solved by using asymmetric bagging [30, 31], a well known and generic method that converts a classifier devoted to well-balanced tasks to unbalanced tasks, which was also discussed in [6] among other approaches. It is worth noticing that we do not solve the issue of false negative, e.g. the fact that among the “no regulation” examples, there might be “regulation” examples that have not been validated yet. The reader interested by this issue is invited to study the works of Cerulo et al. [7] and Mordelet & Vert [6]. Finally, as a fourth contribution, we define and perform three typical numerical studies that can be drawn in order to test a machine learning method devoted to edge prediction: one is a basic test with artificially balanced samples in which we just test the ability of the learning method to obtain good performance; the second one consists of building a regulation predictor in a realistic setting from unbalanced datasets using asymmetric bagging and measuring its ability to discover regulations that were not known before; in the third last study, we proceed in the same way but test the ability of the classifier to label correctly pairs of genes with genes from the training network and genes coming from a new candidate set. In order to assess the performance of the MLN-based approach, we define a pairwise Support Vector Machine (SVM) devoted to ordered pairs of genes and use it as a baseline using a straightforward simplification of the tensor product pairwise kernel. Kernel-based methods as well as first-order logic provide a framework to take into account different sources and features of the data: in this study, two simple definitions of pairwise kernels that combine multiple pairwise kernels expressing heterogeneous information are proposed. While the goal of the study is to take advantage of the heterogeneity of features to describe a pair of genes, we also study the behavior of MLN compared to pairwise SVM in the case of single source of quantitative information such as gene expression.

In order to show the interest of solving these four issues, we have applied our approach to the ID2 genetic regulatory network in human keratinocytes and a new dataset of gene expression using RNA interference. The ID2 protein (Inhibitor of Differentiation 2) acts as a negative regulator of basic helix-loop-helix transcription factors. Previous studies have suggested a potential role for ID2 in epidermis homeostasis reflected by the high expression level of ID2 in proliferating keratinocytes and its down-regulation upon the onset of differentiation [32]. However, the precise implications of ID2 in the process, and in particular its genetic interactions, remain largely unknown. In an attempt to decipher the ID2 genetic regulation network in human keratinocytes, we conducted a transcriptomic analysis by microarray experiments of HaCaT cells presenting stable overexpression or transient knock-down achieved by RNA interference of ID2 expression. As a starting point, we retrieved the regulatory networks associated with the differentially expressed genes in cells with high and low level of ID2 from the Ingenuity Pathway Analysis (IPA) database. We selected a subset of these networks with ontologies of interest for the biologists (cell cycle regulation, cancer, gene expression and signal transduction), merged the corresponding networks and kept only the transcriptional/expression regulations between the genes. The resulting network was finally used to label the couples of genes as a training set.

Methods

Learning directed edges from a set of known regulations

Let$\mathcal{G}$ be the set of genes of interest. We want to learn a function h that takes the descriptors of a gene G₁ and a gene G₂ and predicts if the gene G₁ regulates G₂. Two types of descriptors are considered: descriptors of genes, for instance protein locations within the cell, and relationships between genes reflecting, for instance, if two genes are located on the same chromosome. Let us denote by$\mathcal{X}$ the set of descriptors on genes and by$\mathcal{R}$ the set of relations. A special descriptor expresses the class: given an ordered pair of two genes G₁ and G₂, it is true when G₁ regulates G₂.

In this work we have chosen to use a first-order logic representation, which allows for an easy representation of several objects (here, genes) and their relationships. Facts representing information about objects and their relations are expressed by atomic expressions, called atoms. They are written P(t₁,…,t _n ), where P is a predicate and t₁,…,t _n are terms; a term being either a variable or a constant. In the remainder strings corresponding to constants will start with upper-case letters and strings corresponding to variables with lower-case letters. An atom is said to be ground if all its variables are set to specific values. A ground atom can be true or false, depending of the truth value of the property it expresses. It can therefore be seen as a boolean variable.

Descriptors on genes are thus expressed by expressions of the form Attr(Gname,V), where Attr denotes the attribute, Gname the name of the gene G and V the value taken by G on the descriptor Attr. For instance ProtLoccell(Akt1,Cytoplasm) means that the subcellular localization of AKT1 product protein is the cytoplasm. For sake of simplicity, we have used the name of the gene to define its product. If a gene codes for several proteins, there is no limitation to denote one gene and all its products by different names. A predicate that relates CodesFor(Gname,Pname) is just needed. Relations between genes are expressed by expressions R e l(G₁name,G₂name) where Rel denotes the relation satisfied by genes G₁ and G₂. For instance, Samechro(Cth,Id3) expresses that the genes CTH and ID3 are located on the same chromosome. The property that G₁ regulates G₂ is expressed by the predicate Regulates(G₁name,G₂name). Given two genes G₁ and G₂, we aim to predict whether Regulates(G₁name,G₂name) is true or false. In short, when there is no ambiguity on the genes we write Y = 1 when it is true, and Y = 0 otherwise. We have chosen the probabilistic framework of supervised classification and we search for a classifier h that is based on an estimation of the a posteriori probability P(Y = 1|G₁,G₂). It can be more formally written

$$\begin{array}{l}{h}_{\theta }(x_1,x_2)=\mathit{\text{sgn}}(\widehat{P}(Y=1|{\mathcal{X}}_1=x_1,{\mathcal{X}}_2=x_2,\\ {\mathcal{R}}_{12}=r_{12},\mathcal{B}=b)-\theta ),\end{array}$$

where ${\mathcal{X}}_i=x_i$ represents the description of G _i , ${\mathcal{R}}_{12}=r_{12}$ represents the relations between G₁ and G₂ and $\mathcal{B}=b$ represents the background knowledge. θ is a threshold, whose value will be discussed in the experiments. As shown by this formalization, the learning framework we consider is beyond the classical framework of machine learning in which data is represented only by attributes; it belongs to the ILP (Inductive Logic Programming) domain, a subfield of machine learning that aims at studying relational learning in first-order formalisms [33].

The model we have chosen is a Markov Logic Network, as introduced in [29]. Such a model is defined as a set of weighted first-order formulas. In this paper, we consider only a subset of first-order logic, composed of rules A₁∧…∧A _n  ⇒ Regulates(g₁,g₂), where A₁, …, A _n are atoms. Such restrictions correspond to Horn clauses. The left-hand side of the rule (A₁∧…∧A _n ) is called the body of the rule whereas the right-hand side is called the head of the rule.

Learning a Markov Logic network

Statistical relational learning (SRL) relates to a subfield of machine learning that combines first-order logic rules with probabilistic graphical frameworks. Among the promising approaches to SRL, Markov Logic Networks (MLNs) introduced by Richardson and Domingos [28, 29] are an appealing model. An MLN$\mathcal{M}$ is defined by a set of formulas F = {f _i |i = 1,…,p} and a weight vector w of dimension p, where the clause f _i has an associated weight w _i (negative or positive) that reflects its importance. Therefore, an MLN provides a way of softening first-order logic and encapsulating the weight learning into a probabilistic framework.

A Markov Logic Network together with a finite set of constants C, among which the variables can take their values, defines a Markov Network. This Markov Network can be built by associating a node to each ground atom and by defining a link between two nodes when their corresponding ground atoms occur in the same ground formula. As a consequence, the ground atoms appearing together in a ground clause form a clique in the graph. Let us, for instance, consider the following weighted clause, where g1 and g2 are two variables representing genes:

$$\begin{array}{l}0.3\mathit{\text{Processbio}}(g2,\mathit{\text{Cell\_proliferation}})\\ \wedge \mathit{\text{Processbio}}(g1,\mathit{\text{Negative\_regulation\_of\_cell\_proliferation}})\\ \Rightarrow \mathit{\text{Regulates}}(g1,g2),\end{array}$$

(1)

where the predicate Processbio(Gname,Proc) says that gene G is involved in the biological process annotation Proc of Gene Ontology [34].

Let us suppose that we have two genes A and B. The clause (1) leads to four instantiated clauses, corresponding to the instantiations of g1 and g2 with A and B :

• 1. Processbio(B,Cell_proliferation) ∧ Processbio(A,Negative_regulation_of_cell_proliferation) ⇒ Regulates(A,B)

• 2. Processbio(A,Cell_proliferation) ∧ Processbio(B,Negative_regulation_of_cell_proliferation) ⇒ Regulates(B,A)

• 3. Processbio(A,Cell_proliferation) ∧ Processbio(A,Negative_regulation_of_cell_proliferation) ⇒ Regulates(A,A)

• 4. Processbio(B,Cell_proliferation) ∧ Processbio(B,Negative_regulation_of_cell_proliferation) ⇒ Regulates(B,B)

Variables of the Markov network are the ground atoms occurring in these clauses and they are linked when they occur in the same clause. For instance, the first instantiated clause leads to links between Processbio(B, Cell_proliferation) and Processbio(A,Negative_regulation_of_cell_proliferation), Processbio(B,Cell_proliferation) and Regulates(A,B), and Processbio(A,Negative_regulation_of_ cell_proliferation) and Regulates(A,B). Figure 1 gives the Markov network built from this clause.

Figure 1

Schema of a Markov network built from a MLN with C  = { A,B }. The Markov network has been obtained from the clause (1) by considering a set of two constants A and B. A different color is associated with each instantiated clause.

A world is an assignment of truth values to all possible ground atoms. It is written for short X = x (X denotes the ground atoms and x their truth values). The probability of a world x is given by:

$$P(X=x)=\frac{1}{Z}\mathit{\text{exp}}(\sum _{f_i\in F}{w}_i\times n_i(x)),$$

(2)

where n _i (x) is the number of true groundings of the clause f _i in the world x, and $Z={\sum }_{x}P(X=x)$ is the partition function used for normalization.

For instance, if we consider a world where Processbio(B, Cell_proliferation), Processbio(A,Negative_regulation_ of_cell_proliferation) are true and the other ground atoms are false, then the first instantiated clause is false in this world, whereas all the other instantiated clauses are true (because their premises are false and the logical implication is false). Thus, the number of true groundings of the clause (1) is 3.

For edge prediction, the aim is to infer a classifier for a specific target predicate, given a set of positive and negative examples and background knowledge. We are thus interested in the Conditional Log Likelihood. Given the predicate Y to learn (Y is Regulates), we note examples for this predicate Y _j  = y _j , j = 1,…,n, and Y = y if and only if ∀j,Y _j  = y _j . Given evidence X which corresponds to descriptors of genes, relations between genes and background knowledge, the Conditional Likelihood (CL) can be expressed using the structure of the Markov network:

$$\begin{array}{ll}P(Y=y|X=x)& =\frac{1}{Z_x}exp(\sum _{f_i\in F}{w}_i{n}_i(x,y))\\ =\frac{exp({\sum }_{f_i\in F}{w}_i{n}_i(x,y))}{{\sum }_{z}exp({\sum }_{f_i\in F}{w}_i{n}_i(x,z))},\end{array}$$

(3)

where n _i (x,y) is the number of true groundings of f _i in the world (x,y).

Learning an MLN consists of structure learning, i.e., learning the logical formulas, and parameter learning, i.e. learning the weight of each formula. Completing these two issues simultaneously raises some complexity issues. Therefore, we have chosen to split the learning task into two subtasks. Structure learning can be handled by an inductive logic program (ILP) learner while weight learning can be addressed by maximizing the Conditional Log Likelihood. These subtasks are illustrated in Figure 2.

Figure 2

Schema of the two steps of the MLN learning algorithm.

Learning the candidate rules with Aleph

The system Aleph, developed by Srinivasan [35], is a well known ILP learner that implements the method proposed in [36]. Aleph, like other relational learners, takes as input ground atoms corresponding to positive and negative examples and background knowledge. It also needs language biases, which restrict the set of clauses that can be generated, thus allowing to reduce the size of the search space. These restrictions can correspond to information specified on the predicates, like the place where they occur in the rule, the types of their arguments or the way they will be used (instantiated or not). In our case, we specified that the predicate Regulates occurs in the head of the rule, and the other ones in the body of the rule. Other constraints, such as the maximum number of atoms in a clause or the number of variables, can be defined in order to restrict the form of the rules that can be learned.

The main learning method developed by Aleph, called induce, is sketched in the following:

• 1. Select a positive example not yet covered by a rule

• 2. Build the most specific clause r that covers this example and that satisfies the language biases from the background knowledge. This clause is called the “bottom clause”.

• 3. Search a clause more general than the bottom clause: perform a top-down search (from the most general to the most specific clause) in the search space bounded by r.

• 4. Add the clause with the best score to the current theory and prune redundant clauses.

• 5. Repeat until all positive examples are covered.

Weight learning

Richardson & Domingos [29] proposed performing generative weight learning for a fixed set of clauses by optimizing the pseudo log-likelihood. Several approaches have been proposed for discriminative learning, where the conditional log-likelihood is optimized instead [37, 38]. Huynh & Mooney [39] introduced a weight learning algorithm that targets the case of MLNs containing only non-recursive clauses. In this particular case, each clause contains only one target predicate, thus the grounding of the clauses will contain only one grounded target predicate. This means that the query atoms are all independent given the background atoms. Because of this special assumption on the structure of the model, their approach can perform exact inference when calculating the expected number of true groundings of a clause. Recently, Huynh & Mooney [40] have introduced a discriminative weight learning method based on a max-margin framework.

As we also considered MLNs containing only non-recursive formulas, we used an MAP approach, maximizing the conditional log-likelihood penalized by an ℓ₂ norm constraint:

$$\begin{array}{ll}f(\mathbf{w})& =logP(Y=y|X=x,\mathbf{w})-\lambda \underset{2}{\overset{2}{\left|\right|\mathbf{w}\left|\right|}}\\ =\sum _{j=1}^{n}logP(Y_j=y_j|X=x,\mathbf{w})-\lambda {\left|\right|\mathbf{w}\left|\right|}_2^2,\end{array}$$

with

$$\begin{array}{ll}P(Y_j& =y_j|X=x,\mathbf{w})\\ =\frac{\mathit{\text{exp}}({\sum }_{i\in {\mathcal{F}}_{Y_j}}{w}_i{n}_i(x,y_{[Y_j=y_j]})}{\mathit{\text{exp}}(\sum _{i\in {\mathcal{F}}_{Y_j}}{w}_i{n}_i(x,y_{[Y_j=0]}))+\mathit{\text{exp}}(\sum _{i\in {\mathcal{F}}_{Y_j}}{w}_i{n}_i(x,y_{[Y_j=1]}))},\end{array}$$

where ${\mathcal{F}}_{Y_j}$ is the set of clauses concluding on the target atom Y _j , and $n_i(x,y_{[Y_j=y_j]})$ is the number of true groundings of the ith clause when the atom Y _j is set to the value Y _j . For finding the vector of weights w optimizing this objective function, we used the limited-memory BFGS algorithm [41] implemented in the software ALCHEMY [42].

Materials for inference of the ID2 genetic regulatory network

Data

We conducted a transcriptomic analysis of microarray experiments of HaCaT cells presenting distinct expression levels of ID2. We analyzed three conditions: wild-type cells (wt), stable overexpression (prcID2) or transient knock-down achieved by RNA interference (siID2) of ID2 expression and their corresponding controls. Differentially expressed genes in prcID2 or siID2 versus the corresponding control cells were identified by a t-test analysis using a p-value cut-off of 0.005, a fold-change threshold of 1.5 and Benjamini & Hochberg multiple testing correction [43]. The resulting genes were mapped to genetic networks as defined by IPA tools and the significantly enriched networks associated with cell cycle regulation, cancer, gene expression and signal transduction were merged. In this merged network, only edges and their associated nodes (genes) corresponding to expression/transcriptional regulations were conserved. Genes with incomplete information for all the features were removed. This process led to the selection of a network containing a set of 63 genes, denoted by ${\mathcal{G}}_A$.

In order to use MLNs, we need to describe known properties of genes within the first-order logic setting.

Encoding data

Three low level predicates have been defined to reflect the corresponding experimental conditions:

• The predicate Expwt(Gname,L) states that the expression level of gene G in the wild-type cells is L. In the following results, expression levels values were discretized using equal width discretization [44]: we divided the interval of gene expression values into 5 intervals of equal width.

• Similarly, the predicate that states that the expression level of gene G is L is Expsiid2(Gname,L) when the expression of ID2 has been decreased, and Expprcid2(Gname,L) when it has been increased.

Three other predicates express an increase, a decrease or a lack of change of the expression level between the experience on the wild-type cells and the other experiences: Expmore(Gname,Exp), Expless(Gname,Exp) and Expsame(Gname,Exp), where Exp is either Prcid2 or Siid2.

In order to characterize regulatory interactions, we used other features describing genes. Some of these features concern proteins and not directly genes.

● Physical interaction between proteins: Physical interaction between proteins can provide a hint about the role played by the genes coding for these proteins. In our study, we used the protein interaction data from the IntAct database [45]. We encoded the information of a physical interaction by a predicate containing the name of the genes that are assumed to code the proteins: Inteprot(G₁name,G₂name).

● Subcellular localization of proteins:

● Another interesting information about proteins is their localization in the cell. All proteins were analyzed using the Ingenuity Pathway Analysis Knowledge Base (Ingenuity Systems, www.ingenuity.com), and we encoded the information on the subcellular localization by a predicate ProtLocCell(Gname,Loc) where G is the name of the gene that codes the protein and Loc is the name of the cellular compartment where the protein was found.

● Biological processes:

● We used Gene Ontology [34] to describe the genes by the biological processes in which they are involved. To do so, we have defined a predicate Processbio(Gname,Proc), which says that a gene G is involved in the process Proc.

● Chromosomal location of genes:

● We extracted the genes location on chromosomes and chromosomal bands from the Entrez Gene database [46]. This information is encoded by the predicates Locchro(Gname,Chro) and Locband (Gname,Arm_begin,Band_begin,Arm_end,Band_end). From these predicates, we built two other predicates that we used instead: Samechro(G₁name,G₂name) and Sameband(G₁name,G₂name). These predicates provide information on the proximity between the gene locations of G₁ and G₂.

Choice of a baseline for comparison

In the results, we present a comparison with two pairwise Support Vector Machines (SVMs) used as a baseline approach. Contrary to local classifiers, pairwise classifiers do not need an assumption about known transcription factors: any ordered pair of genes can be processed without any prior. As SVM is built from the definition of a similarity between input data, we need to define a kernel between ordered pairs of genes. We say that two ordered pairs of genes (G₁,G₂) and (G₃,G₄) are similar if the regulator candidate G₁ is similar to the regulator candidate G₃ and similarly, the regulee candidate G₂ is similar to the regulee candidate G₄. This definition requires to first choose a kernel between single data noted k and then writes as:

$$K((G_1,G_2),(G_3,G_4))=k(G_1,G_3)k(G_2,G_4).$$

(4)

This pairwise kernel is the asymmetric version of the kernel proposed in [20, 22] for pairs of proteins to solve supervised protein-protein interaction network inference tasks. Alternative definitions of pairwise kernel have also been proposed, like the metric learning pairwise kernel [47] and the cartesian kernel [48, 49].

For the pairwise kernel defined in (4), when k is chosen to be gaussian and G₁,G₂,G₃,G₄ have a feature vector description, K is also equivalent to a simple unique gaussian kernel built on the concatenation of feature vectors of each pair’s component such as the one proposed in [4]. In the experimental results we present, we defined six gaussian kernels for each feature described previously: gene expressions, differences of gene expression, protein-protein interactions, subcellular localizations, biological processes and chromosomal locations. However the definition proposed in (4) opens the door to different ways of combining the information. We tested two ways of combining kernels that have been proposed in the pairwise SVM framework (see [20] for instance). The first one consists in deriving for each kernel k _i , defined as a kernel between single data, a pairwise kernel k _i and averaging them to build a pairwise kernel noted K _pairwisesum :

$$\begin{array}{ll}{K}_{\mathit{\text{pairwisesum}}}& ((G_1,G_2),(G_3,G_4))\\ =\frac{1}{6}\sum _{i=1}^6{K}_i((G_1,G_2),(G_3,G_4)).\end{array}$$

The second one consists in first averaging the Gaussian kernels and build as final kernel K _sum :

$$K_{\mathit{\text{sum}}}((G_1,G_2),(G_3,G_4))=\bar{k}(G_1,G_3)\bar{k}(G_2,G_4)),$$

where $\bar{k}(G_j,G_k)=\frac{1}{6}{\sum }_{i=1}^6{k}_i(G_j,G_k)$.

Let us notice that kernels are appropriate tools to gather heterogeneous sources of information into the same framework and that combining multiple kernels allows active data integration. Once an SVM is built it is hard to open the “black box” and interpret the decision function.

Results and discussion

Description of the experimental studies

We conducted three numerical studies of the gene regulatory network associated with ID2 in human cells, which are summarized in Table 1.

Table 1 Summary of the three experimental studies

In the first study, we considered the set of 106 regulations provided by Ingenuity in 2007 between the genes in ${\mathcal{G}}_A$, denoted by $R_1^{+}$. All the unknown regulations ($\left|\overline{R_1^{+}}\right|=3863$) were considered as negative examples. The goal of this first study was to test a Markov Logic Network on a well-balanced classification task.

For the second study, we considered the set $R_2^{+}$ of regulations provided by Ingenuity in 2009 for the same set of genes ${\mathcal{G}}_A$. We figured out that 51 new regulations have been discovered by Ingenuity between 2007 and 2009 and we were interested in the prediction task on the updated network. Usual bagging applied to an unbalanced dataset will provide biased classifiers. To build a classifier appropriate for an unbalanced prediction task, we used asymmetric bagging [30, 31].

In supervised classification, asymmetric bagging consists of performing random sampling only on the over-represented class, such that the number of examples in the subsample is equal to the number of examples in the under-represented class. This way, each generated predictor was trained on a balanced dataset. Their predictions on the test set were combined to provide a single prediction. Studies described in [30, 31] have shown that asymmetric bagging provide better results than normal bagging on unbalanced datasets.

In the last study, we solved a network completion task in real conditions. We selected a new set of genes ${\mathcal{G}}_B$ and tried to infer the known regulations between the genes of ${\mathcal{G}}_B$ and ${\mathcal{G}}_A$. Asymmetric bagging was also applied.

The lists of genes in ${\mathcal{G}}_A$ and ${\mathcal{G}}_B$ are given in the Additional file 1 and details on Aleph parameters are available in the Additional file 2. Regarding Alchemy, we used the implementation of the discriminative weights learning procedure and tested different values of the regularization parameter λ.

Evaluation metric

We used area under the ROC (resp. Precision-Recall) curves as evaluation metrics, denoted by AUC-ROC (resp. AUC-PR). These curves were obtained by tuning the threshold θ from 0 to 1 in order to predict regulations from posterior probabilities. It is well known that a ROC curve shows the behavior of the True Positive Rate (also called recall), $\mathit{\text{TPR}}=\frac{\mathit{\text{tp}}}{p}$, according to the value of the False Positive Rate, $\mathit{\text{FPR}}=\frac{\mathit{\text{fp}}}{n}$, while a PR curve assesses the behavior of the precision, $\mathit{\text{Precision}}=\frac{\mathit{\text{tp}}}{\mathit{\text{tp}}+\mathit{\text{fp}}}$, according to the value of the recall. A ROC curve expresses the price to be paid in terms of wrongly predicted negative examples when retrieving correctly a number of positive cases. A PR curve, usually plot in information retrieval tasks, puts emphasis on the confidence of positive predictions. We standardized our precision-recall curves similarly to what was proposed in [50].

Average cross-validation measurements on balanced samples

We first tested the performance of an MLN and compared it to that of a pairwise SVM on a well-balanced classification task. To do that, we subsampled the negative example set and generated subsamples of negative examples of the same size as the positive examples set.

The dataset contains a set $R_1^{+}$ of 106 positive examples of regulations between the genes of ${\mathcal{G}}_A$. We randomly sampled 30 sets of negative examples R1,i−, i = 1,…,30 with $R_{1,i}^{-}\subseteq \overline{R_1^{+}}$, and $\left|R_{1,i}^{-}\right|=\left|R_1^{+}\right|$. Then for each sample we performed a 10-fold cross-validation experiment (10-CV) on $R_1^{+}\cup R_{1,i}^{-}$. In each experiment, we first used Aleph, the ILP tool previously described, as a structure learner on the training set. With Alchemy [42], we learned the weights of the MLN defined by the structure obtained with Aleph and then we performed inference on the test set. For the SVM we used Lib-SVM that we fed with the right definition of the two pairwise kernels computed in Matlab. The bandwidth parameter for each of the six Gaussian kernels was chosen to maximize the entropy of the kernel values. Table 2 shows the averaged AUC-ROC and AUC-PR values obtained within a large range of values of the hyperparameter λ of the MLN while Table 3 contains the same results for bagged SVM for different values of C. As λ, C is a regularization parameter. It controls the importance of the ℓ₁ norm of the slack variables in the dual expression of the loss function. The results of the MLN are not significantly different from those of the best SVM, the “sum” one. These results are very good both in terms of AUC-ROC and AUC-PR. It is also important to notice that neither the MLN nor the SVM are very sensitive to the value of the hyperparameter. However we have noted that λ, the ℓ₂ norm constraint parameter has to be chosen high (larger than 20 to get interesting results).

Table 2 Averaged AUCs for cross-validation measurements on balanced samples using MLNs

Table 3 Averaged AUCs for cross-validation measurements on balanced samples using SVMs

Prediction on the updated graph

In this second study, we addressed a network completion task while keeping the same set of nodes. Two years after the dataset described previously was obtained, the tool Ingenuity was used again to provide an updated set $R_2^{+}$ of regulations between the 63 genes of interest on this date. We noticed that 51 new regulations were discovered by Ingenuity between these two dates. We were therefore interested in the prediction task of the updated graph, i.e. to see if we could retrieve these new regulations from the data of 2007. We used the dataset $R_1^{+}$ from 2007 containing 106 regulations as positive training set and tried to infer the 51 new regulations in $R_2^{+}\setminus R_1^{+}$ using asymmetric bagging. To that end, we randomly sampled 30 negative examples training sets $R_{1,i}^{-}$, i = 1,…,30 with $R_{1,i}^{-}\subseteq \overline{R_1^{+}}\setminus R_2^{+}$ and $\left|R_{1,i}^{-}\right|=\left|R_1^{+}\right|$.

As negative examples correspond here to absences of regulation, the test examples were all positive in this study. We could therefore compute the proportion of regulations which were correctly predicted as positive ones by the classifier with a threshold selected using a validation set in the following way: for each sampling of the negative examples, 2/3 of $R_1^{+}$ and $R_{1,i}^{-}$ were considered for the training set and the remaining regulations were considered for the validation set. We computed each time the F₁-measure obtained on the validation set for different threshold values between 0 and 1:

$$\begin{array}{c}{F}_1=2.\frac{\mathit{\text{Precision}}\ast \mathit{\text{Recall}}}{\mathit{\text{Precision}}+\mathit{\text{Recall}}}.\end{array}$$

We selected the threshold maximizing the averaged F₁-measure, that is the value maximizing precision and recall at the same time.

Then, for each sampling, we applied the predictor learned on the training set to the 51 new regulations. We averaged the predictions obtained and used the selected threshold value to compute the true positive rate (TPR). The TPR values obtained using bagged MLN and bagged pairwise SVMs respectively are given in Table 4. Again results using the MLN are very good, showing that it is possible to predict new regulations from an existing corpus at a given time. The performance of pairwise SVMs are inverted, the “pairwise sum” SVM achieving the same results as the MLN.

Table 4 Prediction of regulations on the updated network

Prediction with a new set of genes

For the third statistical analysis, we addressed a network completion task when new candidate nodes are added. We used a dataset refined in the biology laboratory. 209 high confidence differentially expressed genes in prcID2 versus the corresponding control cells were identified. From these genes, we selected 37 genes that were not part of ${\mathcal{G}}_A$ and for which we had an annotation for each predicate. These genes were also chosen from the ones that had at least one regulation link with one of the genes from ${\mathcal{G}}_A$ or with one gene of this new set. From these 37 genes, we selected a subset of 24 genes, called ${\mathcal{G}}_B$, that had at least a biological process annotation from GO in common with genes from ${\mathcal{G}}_A$. The goal of this study was to try to complete a known network using an additional set of candidates genes, which is usually the problem of interest for the biologists. We used Ingenuity to retrieve the known regulations between genes from ${\mathcal{G}}_A$ and ${\mathcal{G}}_B$, being aware that when no regulation is mentioned in the literature, it does not mean that it does not exist but only that it has not been discovered yet. We called this set $R_3^{+}$.

We used the set $R_2^{+}$ of 157 regulations from the dataset ${\mathcal{G}}_A$ as the positive examples training set and used asymmetric bagging with 30 samples on R₂. For each sampling we applied the predictor on the sets $R_3^{+}$ and $\overline{R_3^{+}}$, using descriptors on both set of genes ${\mathcal{G}}_A$ and ${\mathcal{G}}_B$. We obtained score predictions for each interaction between one gene of ${\mathcal{G}}_A$ and one gene of ${\mathcal{G}}_B$. Table 5 reports the AUC values computed for bagged MLNs and bagged SVMs. The ROC and PR curves obtained for the hyper-parameters associated with the best values of AUC-ROC are represented in the Additional file 3.

Table 5 Prediction of regulations between the set of genes ${\mathcal{G}}_{\mathit{A}}$ and ${\mathcal{G}}_{\mathit{B}}$

Although each predictor was trained on a balanced dataset, with same numbers of positive and negative examples of regulation, this test was made under real conditions: we considered the whole set of positive ($\left|R_3^{+}\right|=55$) and negative examples ($\left|\overline{R_3^{+}}\right|=2969$) to assess the performance in prediction. On the test-training interactions, the predictor with bagged MLNs performed quite well, showing an AUC-ROC of about 0.73. This was really a very good result which implies low degradation in performance especially for the false positive rate that only slightly increases. The AUC values obtained with bagged MLNs are above the values obtained with the two bagged SVMs. We performed a statistical test in order to compare the AUC-ROC values obtained with the different classifiers. We used the non-parametric test on Mann Whitney statistics developed by [51] and the implementation provided by the R package pROC [52]. The obtained p-values are given in the Additional file 4. We observe from this results that the p-values are less than 0.05 and therefore that the AUC-ROC values of bagged MLNs and bagged pairwise sum SVMs are significantly different. Regarding the comparison between bagged MLNs and bagged sum SVMs, the difference between AUC-ROC values is not significant, indicating similar predictive performance.

AUC-PR of bagged MLNs outperforms the best pairwise SVM. Therefore in a real prediction task, e.g. a network completion task, MLN exhibits a very interesting behaviour, even if the AUC-PR still needs to be increased.

In Table 6, we also present the results obtained for this task when using only gene expression levels as input descriptors. First we notice that, compared to the results in Table 5, the performance of both approaches, SVMs and MLNs, diminish, showing that the additional descriptors play an important role. Second, in terms of AUC-ROC, the performance of bagged pairwise SVMs and bagged MLNs are similar. The p-values obtained by performing a statistical test for the comparison of these AUC-ROC values are greater than 0.3 for all hyper-parameter values. We can therefore deduce that the difference of AUC-ROC between bagged MLNs and bagged SVMs in terms of AUC-ROC is not significant. Third, for all the approaches, AUC-PR values are very low, but bagged SVMs now outperform bagged MLNs. As expected, SVMs are more suitable for dealing with numerical data. We therefore recommend to use the MLN method when different sources of gene descriptions such as symbolic and discrete descriptions are available.

Table 6 Prediction of regulations between the set of genes ${\mathcal{G}}_{\mathit{A}}$ and ${\mathcal{G}}_{\mathit{B}}$ when using only gene expression data as descriptors

To conclude, we have shown in this section that bagged sum SVM performs well in Task1 and Task3, while bagged pairwise sum SVM performs well in Task2. Contrary to the SVM classifiers, MLNs behaved well in the three tasks. Now another interesting criterion to choose a method for network inference is to measure its ability to provide insights on the taken decisions.

Resulting logical rules

In addition to the capacity of the built classifier to suggest new regulations, MLNs present the advantage of providing a set of weighted rules that the biologist can check. In general, Aleph learned between 30 and 50 rules for each run, these rules being composed of up to five predicates. This comes from the choice of the parameters of Aleph as described in the Additional file 2. The analysis of the rules that have obtained a high weight shows that some of them exhibit relevant patterns. Among the rules inferred by Aleph, here is an example of four rules which were associated to a high weight in the numerical tests:

• 1. ProtLoccell(g ₂,Plasma_membrane) ∧ Expsiid2(g ₂,Level3) ∧ Expsiid2(g ₁,Level3) ⇒ Regulates(g ₁,g ₂)

• 2. Processbio(g ₂,Cell_proliferation) ∧ Processbio(g ₁,Negative_regulation_of_cel_proliferation) ⇒ Regulates(g ₁,g ₂)

• 3. Expsiid2(g ₁,Level3) ∧ Expprcid2(g ₁,Level4) ∧ Expsiid2(g ₂,Level4) ∧ Expprcid2(g ₂,Level5) ⇒ Regulates(g ₁,g ₂)

• 4. Expprcid2(g ₁,Level5) ∧ Expwt(g ₂,Level2) ∧ Expprcid2(g ₂,Level4) ⇒ Regulates(g ₁,g ₂)

The first rule means that a gene overexpressed in transient knock-down of ID2 regulates overexpressed genes in the same condition and that code for proteins in plasma membrane. Obviously, this rule alone is far too general but within a set of rules with positive and negative weights, it brings a piece of evidence for regulation. The second rule may seem trivial but it has been retrieved from data: it says that genes involved in negative regulation of cell proliferation regulate genes involved in cell proliferation. The next rule means that an increase of the expressions of G₁ and G₂ in the condition of over expression of ID2 compared to transient knock-down of ID2 indicates a regulation between G₁ and G₂. Regarding the last rule, it indicates that a high expression value of G₁ in the prcID2 condition and the increase of the expression of G₂ between wild-type condition and prcID2 implies the existence of a regulation between these two genes.

These rules are examples of what has been obtained in a first attempt to build a whole strategy to get a supervised edge predictor. However the quality of the learnt rules strongly depends on the nature of the chosen predicates and the ILP learning phase. We notice that a substantial improvement can be reached in terms of rules if the biologist makes explicit some constraints on the rules. For instance, one might want rules that include at least relations on both input genes in their premises. We will favor this research direction in the future.

Another information that can be extracted from the learnt MLN concerns the statistics of presence of some of the predicates in the premises of the rules. In our experimental studies, chromosomal location of genes did not appear as an important property to conclude about regulation.

Conclusions

Recent years have witnessed the preeminence of numerical and statistical machine learning tools in computational biology. Among them, kernel-based methods present numerous advantages, including the ability to deal with heterogeneous sources of information by encoding them into similarities (kernels). On top of that, multiple kernel learning allows to select sources of informations thought the learning of sparse linear combination of kernels [19, 53, 54]. However kernel-based methods remain black boxes: using non linear kernels, the decision function learnt with a SVM is not at all interpretable. This is an inherent drawback of SVMs because biologists are generally not only interested in the prediction made by a classifier but also in the reason why such an example has been labeled in a given way.

This work explores another direction through a new hybrid tool based on first-order logic representation and probabilistic graphical modeling. Once learnt, a MLN provides a weighted set of rules that conclude on the target predicate, here the regulates predicate. To our knowledge this work is the first application of MLN to gene regulatory network inference and one of the very first real applications of MLN on noisy and medium scale biological data. As described in the previous sections, learning a MLN involves several steps including data encoding, choice of constraints and hyper-parameters in the ILP learner and the weight learner as well as an appropriate learning protocol scheme for achieving the learning task. All these steps require a high level of collaboration between biologists and computer scientists which is facilitated by the common language of first-order logic. Therefore, in one hand, the encoding process can be seen as a limitation since each new application requires specific work about the choice and the definition of the predicates to be used. Compared to the kernel design, this step is expensive. However, on the other hand, it produces a corpus of interpretable facts and rules encoding the nature of the relationship between genes that the biologist can inspect. Moreover, it is worth pointing out the fact that it is relatively easy in this context to impose known rules or to perform incremental learning at the level of the rule learner. There is also a lot of relevant information that can be made available that we did not incorporate to describe genes. For instance, adding knowledge of regulatory motifs of genes and DNA-binding sites of regulatory proteins, could improve the performance of the predictor. This means that a proper representation of sequences should be described either directly in first-order logic as it was done in [55], or using an extension of first-order logic to sequence variables like those of [56]. This is certainly a direction to be explored in future works.

Another issue is scalability to larger networks composed of thousands of genes. This would be a concern for pairwise kernel-based methods for instance for the later task to compute the Gram matrix between training and test data. For MLN, scaling to a larger number of genes like thousands of genes should be made possible using the latest improvement in MLN learning implemented in FELIX [57] using dual composition.

Another interesting question is to compare decision trees with MLNs. Decision trees are usually built from attribute-value representations but have been extended to first-order logic in [58]. They also provide a set of interpretable rules but in a less general form than in MLNs. In a decision tree, rules are factorisable and a given example to be classified will only satisfy one rule. On the contrary, a MLN devoted to supervised classification a given example can satisfy many rules. Interestingly, combining decision trees to learn compact representations of MLNs has been recently proposed in [59].

Finally the biologist interested in the ID2 genetic regulatory network in human keratinocytes gets two main results from this work additionally to a set of facts and rules describing the network. First, learning such a supervised pairwise classifier can be seen as a cross-validation of both experiments and existing literature. The ability of the learning algorithm to build a good edge prediction tool shows indeed that text-mining and careful curation can produce networks that are consistent. Inversely, the experimental data measured in the wet laboratory are proven to make sense. Second, the last in silico study can provide a list of predicted regulations with new candidate genes, some of them being known but some of the others, considered currently as false positive, may involve new regulators and new targets. This calls for an experimental wet lab validation to test the relevance of the potential new regulations.