Learning a Markov Logic network for supervised gene regulatory network inference

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.

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.

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

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.

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.

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.

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.

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.

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

Choice of a baseline for comparison

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].

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.

Description of the 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 λ.

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.

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|$.

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

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^{+}$.

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.

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

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.