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  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  <cit> .
three main families of inference algorithms have been developed so far:  unsupervised model-free approaches that use information theory to extract a non-oriented graph of dependence between variables,  unsupervised reverse-modeling approaches that model the network behavior as a  system  <cit>  and  supervised edge prediction approaches that focus on the graph of regulation and only predict the presence/absence of regulations  <cit> .
in the first family, relevance networks like aracne  <cit> , clr  <cit>  and td-aracne  <cit>  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   <cit>  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  <cit>  has also proven to outperform competitors in inferring gene regulatory networks in recent dream competitions.
bayesian networks  <cit>  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  <cit> .
taking a different angle, supervised edge prediction methods build a decision function that associates a class label to a pair of vertices  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  <cit> .
the principle underlying  <cit>  is to build pairwise support vector machines  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  to deal with multiple sources of information.
in  <cit> , another point of view is taken: local models  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  <cit>  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.
<cit>  using gene expression as unique descriptor and further developed by mordelet et al.
with the sirene method  <cit> .
similarly to  <cit> , 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  <cit> ; others adjust the probability of being positive estimated by a classifier trained on positive and unlabeled examples  <cit> .
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  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 id <dig> 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  <cit> , we choose a markov logic network   <cit>  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  <cit> , a well known and generic method that converts a classifier devoted to well-balanced tasks to unbalanced tasks, which was also discussed in  <cit>  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.
<cit>  and mordelet & vert  <cit> .
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  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 id <dig> genetic regulatory network in human keratinocytes and a new dataset of gene expression using rna interference.
the id <dig> protein  acts as a negative regulator of basic helix-loop-helix transcription factors.
previous studies have suggested a potential role for id <dig> in epidermis homeostasis reflected by the high expression level of id <dig> in proliferating keratinocytes and its down-regulation upon the onset of differentiation  <cit> .
however, the precise implications of id <dig> in the process, and in particular its genetic interactions, remain largely unknown.
in an attempt to decipher the id <dig> 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 id <dig> expression.
as a starting point, we retrieved the regulatory networks associated with the differentially expressed genes in cells with high and low level of id <dig> from the ingenuity pathway analysis  database.
we selected a subset of these networks with ontologies of interest for the biologists , 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 g be the set of genes of interest.
we want to learn a function h that takes the descriptors of a gene g <dig> and a gene g <dig> and predicts if the gene g <dig> regulates g <dig>  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 x the set of descriptors on genes and by r the set of relations.
a special descriptor expresses the class: given an ordered pair of two genes g <dig> and g <dig>  it is true when g <dig> regulates g <dig>
in this work we have chosen to use a first-order logic representation, which allows for an easy representation of several objects  and their relationships.
facts representing information about objects and their relations are expressed by atomic expressions, called atoms.
they are written p, where p is a predicate and t <dig> …,tn 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, 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 means that the subcellular localization of akt <dig> 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 is just needed.
relations between genes are expressed by expressions rel where rel denotes the relation satisfied by genes g <dig> and g <dig>  for instance, samechro expresses that the genes cth and id <dig> are located on the same chromosome.
the property that g <dig> regulates g <dig> is expressed by the predicate regulates.
given two genes g <dig> and g <dig>  we aim to predict whether regulates is true or false.
in short, when there is no ambiguity on the genes we write y =  <dig> when it is true, and y =  <dig> 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. it can be more formally written
hθ=sgn−θ),
where xi=xi represents the description of gi, r12=r <dig> represents the relations between g <dig> and g <dig> and 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  domain, a subfield of machine learning that aims at studying relational learning in first-order formalisms  <cit> .
the model we have chosen is a markov logic network, as introduced in  <cit> .
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 a1∧…∧an ⇒ regulates, where a <dig>  …, an are atoms.
such restrictions correspond to horn clauses.
the left-hand side of the rule  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  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  introduced by richardson and domingos  <cit>  are an appealing model.
an mln m is defined by a set of formulas f = {fi|i =  <dig> …,p} and a weight vector w of dimension p, where the clause fi has an associated weight wi  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 g <dig> and g <dig> are two variables representing genes:
<dig> processbio∧processbio⇒regulates,
where the predicate processbio says that gene g is involved in the biological process annotation proc of gene ontology  <cit> .
let us suppose that we have two genes a and b. the clause  leads to four instantiated clauses, corresponding to the instantiations of g <dig> and g <dig> with a and b :
<dig>  processbio ∧ processbio ⇒ regulates
<dig>  processbio ∧ processbio ⇒ regulates
<dig>  processbio ∧ processbio ⇒ regulates
<dig>  processbio ∧ processbio ⇒ regulates
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 and processbio, processbio and regulates, and processbio and regulates.
figure  <dig> gives the markov network built from this clause.
a world is an assignment of truth values to all possible ground atoms.
it is written for short x = x .
the probability of a world x is given by:
p=1zexp),
where ni is the number of true groundings of the clause fi in the world x, and z=∑xp is the partition function used for normalization.
for instance, if we consider a world where processbio, processbio 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 .
thus, the number of true groundings of the clause  is  <dig>
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 , we note examples for this predicate yj = yj, j =  <dig> …,n, and y = y if and only if ∀j,yj = yj.
given evidence x which corresponds to descriptors of genes, relations between genes and background knowledge, the conditional likelihood  can be expressed using the structure of the markov network:
p=1zxexp)=exp)∑zexp),
where ni is the number of true groundings of fi in the world .
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  learner while weight learning can be addressed by maximizing the conditional log likelihood.
these subtasks are illustrated in figure  <dig>
learning the candidate rules with aleph
the system aleph, developed by srinivasan  <cit> , is a well known ilp learner that implements the method proposed in  <cit> .
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 .
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:
<dig>  select a positive example not yet covered by a rule
<dig>  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”.
<dig>  search a clause more general than the bottom clause: perform a top-down search  in the search space bounded by r.
<dig>  add the clause with the best score to the current theory and prune redundant clauses.
<dig>  repeat until all positive examples are covered.
weight learning
richardson & domingos  <cit>  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  <cit> .
huynh & mooney  <cit>  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  <cit>  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 ℓ <dig> norm constraint:
f=logp−λ||w||22=∑j=1nlogp−λ||w|| <dig>
with
p=expexp)+exp),
where fyj is the set of clauses concluding on the target atom yj, and ni is the number of true groundings of the ith clause when the atom yj is set to the value yj.
for finding the vector of weights w optimizing this objective function, we used the limited-memory bfgs algorithm  <cit>  implemented in the software alchemy  <cit> .
materials for inference of the id <dig> genetic regulatory network
data
we conducted a transcriptomic analysis of microarray experiments of hacat cells presenting distinct expression levels of id <dig>  we analyzed three conditions: wild-type cells , stable overexpression  or transient knock-down achieved by rna interference  of id <dig> expression and their corresponding controls.
differentially expressed genes in prcid <dig> or siid <dig> versus the corresponding control cells were identified by a t-test analysis using a p-value cut-off of  <dig> , a fold-change threshold of  <dig>  and benjamini & hochberg multiple testing correction  <cit> .
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  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  <dig> genes, denoted by ga.
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 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  <cit> : we divided the interval of gene expression values into  <dig> intervals of equal width.
● similarly, the predicate that states that the expression level of gene g is l is expsiid <dig> when the expression of id <dig> has been decreased, and expprcid <dig> 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, expless and expsame, where exp is either prcid <dig> or siid <dig>
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  <cit> .
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.
● 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 , and we encoded the information on the subcellular localization by a predicate protloccell 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  <cit>  to describe the genes by the biological processes in which they are involved.
to do so, we have defined a predicate processbio, 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  <cit> .
this information is encoded by the predicates locchro and locband .
from these predicates, we built two other predicates that we used instead: samechro and sameband.
these predicates provide information on the proximity between the gene locations of g <dig> and g <dig>
choice of a baseline for comparison
in the results, we present a comparison with two pairwise support vector machines  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  and  are similar if the regulator candidate g <dig> is similar to the regulator candidate g <dig> and similarly, the regulee candidate g <dig> is similar to the regulee candidate g <dig>  this definition requires to first choose a kernel between single data noted k and then writes as:
k,)=kk.
this pairwise kernel is the asymmetric version of the kernel proposed in  <cit>  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  <cit>  and the cartesian kernel  <cit> .
for the pairwise kernel defined in , when k is chosen to be gaussian and g <dig> g <dig> g <dig> g <dig> 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  <cit> .
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  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 .
the first one consists in deriving for each kernel ki, defined as a kernel between single data, a pairwise kernel ki and averaging them to build a pairwise kernel noted kpairwisesum:
kpairwisesum,)=16∑i=16ki,).
the second one consists in first averaging the gaussian kernels and build as final kernel ksum:
ksum,)=k¯k¯),
where k¯=16∑i=16ki.
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
description of the experimental studies
we conducted three numerical studies of the gene regulatory network associated with id <dig> in human cells, which are summarized in table  <dig>
we conducted three experimental studies on the gene regulatory network associated with id <dig> in human cells.
in the table, 10-cv means cross-validation  <dig> times and ab means asymmetric bagging.
in the first study, we considered the set of  <dig> regulations provided by ingenuity in  <dig> between the genes in ga, denoted by r1+.
all the unknown regulations  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 r2+ of regulations provided by ingenuity in  <dig> for the same set of genes ga. we figured out that  <dig> new regulations have been discovered by ingenuity between  <dig> and  <dig> 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  <cit> .
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  <cit>  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 gb and tried to infer the known regulations between the genes of gb and ga. asymmetric bagging was also applied.
the lists of genes in ga and gb are given in the additional file  <dig> and details on aleph parameters are available in the additional file  <dig>  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  curves as evaluation metrics, denoted by auc-roc .
these curves were obtained by tuning the threshold θ from  <dig> to  <dig> in order to predict regulations from posterior probabilities.
it is well known that a roc curve shows the behavior of the true positive rate , tpr=tpp, according to the value of the false positive rate, fpr=fpn, while a pr curve assesses the behavior of the precision, precision=tptp+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  <cit> .
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 r1+ of  <dig> positive examples of regulations between the genes of ga. we randomly sampled  <dig> sets of negative examples r <dig> i−, i =  <dig> …, <dig> with r <dig> i−⊆r1+¯, and |r <dig> i−|=|r1+|.
then for each sample we performed a 10-fold cross-validation experiment  on r1+∪r <dig> i−.
in each experiment, we first used aleph, the ilp tool previously described, as a structure learner on the training set.
with alchemy  <cit> , 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  <dig> shows the averaged auc-roc and auc-pr values obtained within a large range of values of the hyperparameter λ of the mln while table  <dig> contains the same results for bagged svm for different values of c. as λ, c is a regularization parameter.
it controls the importance of the ℓ <dig> 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 ℓ <dig> norm constraint parameter has to be chosen high .
the table reports the averaged auc values and standard deviations obtained with mlns for thirty ten folds cross-validation experiments conducted on a regulatory network between the genes in ga. the results are reported for different values of the regularization parameter λ.
the table reports the averaged auc values and standard deviations obtained with svms for thirty ten folds cross-validation experiments conducted on a regulatory network between the genes in ga. the column “sum” shows the results when the pairwise kernel is derived from the sum of genomic kernels, while the column “pairwise sum” shows the results obtained using the sum of pairwise kernels derived from each genomic kernel.
the results are reported for different values of the regularization parameter c.
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 r2+ of regulations between the  <dig> genes of interest on this date.
we noticed that  <dig> 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  <dig>  we used the dataset r1+ from  <dig> containing  <dig> regulations as positive training set and tried to infer the  <dig> new regulations in r2+∖r1+ using asymmetric bagging.
to that end, we randomly sampled  <dig> negative examples training sets r <dig> i−, i =  <dig> …, <dig> with r <dig> i−⊆r1+¯∖r2+ and |r <dig> i−|=|r1+|.
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/ <dig> of r1+ and r <dig> i− were considered for the training set and the remaining regulations were considered for the validation set.
we computed each time the f1-measure obtained on the validation set for different threshold values between  <dig> and 1:
() f1= <dig> precision∗recallprecision+recall.
we selected the threshold maximizing the averaged f1-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  <dig> new regulations.
we averaged the predictions obtained and used the selected threshold value to compute the true positive rate .
the tpr values obtained using bagged mln and bagged pairwise svms respectively are given in table  <dig>  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.
this table lists the true positive rates  obtained for the prediction of regulations in r2+∖r1+ from r <dig> using asymmetric bagging with  <dig> samples for bagged mlns and bagged pairwise svms.
the tpr values were obtained using a threshold maximizing the averaged f1-measure on a validation set.
notations are given in tables  <dig> and  <dig>
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.
<dig> high confidence differentially expressed genes in prcid <dig> versus the corresponding control cells were identified.
from these genes, we selected  <dig> genes that were not part of ga 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 ga or with one gene of this new set.
from these  <dig> genes, we selected a subset of  <dig> genes, called gb, that had at least a biological process annotation from go in common with genes from ga. 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 ga and gb, 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 r3+.
we used the set r2+ of  <dig> regulations from the dataset ga as the positive examples training set and used asymmetric bagging with  <dig> samples on r <dig>  for each sampling we applied the predictor on the sets r3+ and r3+¯, using descriptors on both set of genes ga and gb.
we obtained score predictions for each interaction between one gene of ga and one gene of gb.
table  <dig> 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  <dig>
this table lists the auc-roc and auc-pr values obtained for the prediction of regulations betweenga and gb for bagged mlns and bagged pairwise svms, with notations given in tables  <dig> and  <dig>  these results were obtained using asymmetric bagging with  <dig> samples on the set r <dig>
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  and negative examples  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  <dig> .
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  <cit>  and the implementation provided by the r package proc  <cit> .
the obtained p-values are given in the additional file  <dig>  we observe from this results that the p-values are less than  <dig>  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  <dig>  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  <dig>  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  <dig>  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.
this table presents the results obtained with bagged mlns and bagged pairwise svms on the third task when using only gene expression data as gene descriptors.
to conclude, we have shown in this section that bagged sum svm performs well in task <dig> and task <dig>  while bagged pairwise sum svm performs well in task <dig>  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  <dig> and  <dig> 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  <dig>  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:
<dig>  protloccell ∧ expsiid2 ∧ expsiid2 ⇒ regulates
<dig>  processbio ∧ processbio ⇒ regulates
<dig>  expsiid2 ∧ expprcid2 ∧ expsiid2 ∧ expprcid2 ⇒ regulates
<dig>  expprcid2 ∧ expwt ∧ expprcid2 ⇒ regulates
the first rule means that a gene overexpressed in transient knock-down of id <dig> 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 <dig> and g <dig> in the condition of over expression of id <dig> compared to transient knock-down of id <dig> indicates a regulation between g <dig> and g <dig>  regarding the last rule, it indicates that a high expression value of g <dig> in the prcid <dig> condition and the increase of the expression of g <dig> between wild-type condition and prcid <dig> 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 .
on top of that, multiple kernel learning allows to select sources of informations thought the learning of sparse linear combination of kernels  <cit> .
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  <cit> , or using an extension of first-order logic to sequence variables like those of  <cit> .
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  <cit>  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  <cit> .
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  <cit> .
finally the biologist interested in the id <dig> 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.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
cb checked the encoding of the data in first-order logic, built the new datasets, developed codes for the whole pipeline and tested all the methods.
fab and cv conceived the methodological part of the project, mad conceived its biological part about the id <dig> network.
dc and mad did the experiments in the wet laboratory and provide a preliminary analysis of the data.
jd encoded the data and the background knowledge into first-order logic.
cv supervised the use of ilp for rule learning and fab supervised the mln approach, the asymmetric bagging part and designed with cb the baseline svm.
mad participated to design the performance assessment with cb and fab.
cb and fab drafted the document with the help of cv and mad.
all authors read and approved the final manuscript.
supplementary material
additional file 1
list of gene symbols used in the study.
click here for file
additional file 2
choice of aleph parameters.
click here for file
additional file 3
roc and pr curves obtained for the prediction of regulations between setsga and gb.
the figure represents the roc and pr curves obtained with bagged mlns, bagged svms  and bagged svms  using the hyperparameters associated with the best auc-roc values.
click here for file
additional file 4
p-values of the non-parametric test based on mann whitney statistics to compare auc-roc obtained with bagged mlns and bagged pairwise svms.
click here for file
acknowledgements
we thank xavier gidrol  and vincent frouin  for their fruitful comments at the beginning of the study.
this work was supported by the agence nationale de la recherche .
the work of cb and fab was completed using funding from .

 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 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. three main families of inference algorithms have been developed so far: unsupervised model-free approaches that use information theory to extract a non-oriented graph of dependence between variables, unsupervised reverse-modeling approaches that model the network behavior as a system and supervised edge prediction approaches that focus on the graph of regulation and only predict the presence/absence of regulations. in the first family, relevance networks like aracne, clr and td-aracne 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 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 has also proven to outperform competitors in inferring gene regulatory networks in recent dream competitions. bayesian networks 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. taking a different angle, supervised edge prediction methods build a decision function that associates a class label to a pair of vertices 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. the principle underlying is to build pairwise support vector machines 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 to deal with multiple sources of information. in, another point of view is taken: local models 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 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. using gene expression as unique descriptor and further developed by mordelet et al. with the sirene method. similarly to, 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; others adjust the probability of being positive estimated by a classifier trained on positive and unlabeled examples. 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 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 id 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, we choose a markov logic network 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, a well known and generic method that converts a classifier devoted to well-balanced tasks to unbalanced tasks, which was also discussed in 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. and mordelet & vert. 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 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 id genetic regulatory network in human keratinocytes and a new dataset of gene expression using rna interference. the id protein acts as a negative regulator of basic helix-loop-helix transcription factors. previous studies have suggested a potential role for id in epidermis homeostasis reflected by the high expression level of id in proliferating keratinocytes and its down-regulation upon the onset of differentiation. however, the precise implications of id in the process, and in particular its genetic interactions, remain largely unknown. in an attempt to decipher the id 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 id expression. as a starting point, we retrieved the regulatory networks associated with the differentially expressed genes in cells with high and low level of id from the ingenuity pathway analysis database. we selected a subset of these networks with ontologies of interest for the biologists, 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 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 x the set of descriptors on genes and by 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 and their relationships. facts representing information about objects and their relations are expressed by atomic expressions, called atoms. they are written p, where p is a predicate and t …,tn 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, 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 means that the subcellular localization of akt 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 is just needed. relations between genes are expressed by expressions rel where rel denotes the relation satisfied by genes g and g for instance, samechro expresses that the genes cth and id are located on the same chromosome. the property that g regulates g is expressed by the predicate regulates. given two genes g and g we aim to predict whether regulates is true or false. in short, when there is no ambiguity on the genes we write y =  when it is true, and y =  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. it can be more formally written hθ=sgn−θ), where xi=xi represents the description of gi, r12=r represents the relations between g and g and 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 domain, a subfield of machine learning that aims at studying relational learning in first-order formalisms. the model we have chosen is a markov logic network, as introduced in. 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 a1∧…∧an ⇒ regulates, where a …, an are atoms. such restrictions correspond to horn clauses. the left-hand side of the rule 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 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 introduced by richardson and domingos are an appealing model. an mln m is defined by a set of formulas f = {fi|i =  …,p} and a weight vector w of dimension p, where the clause fi has an associated weight wi 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 g and g are two variables representing genes: processbio∧processbio⇒regulates, where the predicate processbio says that gene g is involved in the biological process annotation proc of gene ontology. let us suppose that we have two genes a and b. the clause leads to four instantiated clauses, corresponding to the instantiations of g and g with a and b: processbio ∧ processbio ⇒ regulates processbio ∧ processbio ⇒ regulates processbio ∧ processbio ⇒ regulates processbio ∧ processbio ⇒ regulates 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 and processbio, processbio and regulates, and processbio and regulates. figure gives the markov network built from this clause. a world is an assignment of truth values to all possible ground atoms. it is written for short x = x. the probability of a world x is given by: p=1zexp), where ni is the number of true groundings of the clause fi in the world x, and z=∑xp is the partition function used for normalization. for instance, if we consider a world where processbio, processbio 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. thus, the number of true groundings of the clause is 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, we note examples for this predicate yj = yj, j =  …,n, and y = y if and only if ∀j,yj = yj. given evidence x which corresponds to descriptors of genes, relations between genes and background knowledge, the conditional likelihood can be expressed using the structure of the markov network: p=1zxexp)=exp)∑zexp), where ni is the number of true groundings of fi in the world. 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 learner while weight learning can be addressed by maximizing the conditional log likelihood. these subtasks are illustrated in figure learning the candidate rules with aleph the system aleph, developed by srinivasan, is a well known ilp learner that implements the method proposed in. 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. 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: select a positive example not yet covered by a rule 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”. search a clause more general than the bottom clause: perform a top-down search in the search space bounded by r. add the clause with the best score to the current theory and prune redundant clauses. repeat until all positive examples are covered. weight learning richardson & domingos 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. huynh & mooney 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 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: f=logp−λ||w||22=∑j=1nlogp−λ||w|| with p=expexp)+exp), where fyj is the set of clauses concluding on the target atom yj, and ni is the number of true groundings of the ith clause when the atom yj is set to the value yj. for finding the vector of weights w optimizing this objective function, we used the limited-memory bfgs algorithm implemented in the software alchemy. materials for inference of the id genetic regulatory network data we conducted a transcriptomic analysis of microarray experiments of hacat cells presenting distinct expression levels of id we analyzed three conditions: wild-type cells, stable overexpression or transient knock-down achieved by rna interference of id expression and their corresponding controls. differentially expressed genes in prcid or siid versus the corresponding control cells were identified by a t-test analysis using a p-value cut-off of, a fold-change threshold of and benjamini & hochberg multiple testing correction. 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 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 genes, denoted by ga. 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 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: we divided the interval of gene expression values into intervals of equal width. ● similarly, the predicate that states that the expression level of gene g is l is expsiid when the expression of id has been decreased, and expprcid 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, expless and expsame, where exp is either prcid or siid 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. 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. ● 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, and we encoded the information on the subcellular localization by a predicate protloccell 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 to describe the genes by the biological processes in which they are involved. to do so, we have defined a predicate processbio, 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. this information is encoded by the predicates locchro and locband. from these predicates, we built two other predicates that we used instead: samechro and sameband. 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 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 and 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,)=kk. this pairwise kernel is the asymmetric version of the kernel proposed in 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 and the cartesian kernel. for the pairwise kernel defined in, 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. 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 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. the first one consists in deriving for each kernel ki, defined as a kernel between single data, a pairwise kernel ki and averaging them to build a pairwise kernel noted kpairwisesum: kpairwisesum,)=16∑i=16ki,). the second one consists in first averaging the gaussian kernels and build as final kernel ksum: ksum,)=k¯k¯), where k¯=16∑i=16ki. 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 we conducted three numerical studies of the gene regulatory network associated with id in human cells, which are summarized in table we conducted three experimental studies on the gene regulatory network associated with id in human cells. in the table, 10-cv means cross-validation times and ab means asymmetric bagging. in the first study, we considered the set of regulations provided by ingenuity in between the genes in ga, denoted by r1+. all the unknown regulations 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 r2+ of regulations provided by ingenuity in for the same set of genes ga. we figured out that new regulations have been discovered by ingenuity between and 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. 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 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 gb and tried to infer the known regulations between the genes of gb and ga. asymmetric bagging was also applied. the lists of genes in ga and gb are given in the additional file and details on aleph parameters are available in the additional file 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 curves as evaluation metrics, denoted by auc-roc. these curves were obtained by tuning the threshold θ from to in order to predict regulations from posterior probabilities. it is well known that a roc curve shows the behavior of the true positive rate, tpr=tpp, according to the value of the false positive rate, fpr=fpn, while a pr curve assesses the behavior of the precision, precision=tptp+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. 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 r1+ of positive examples of regulations between the genes of ga. we randomly sampled sets of negative examples r i−, i =  …, with r i−⊆r1+¯, and |r i−|=|r1+|. then for each sample we performed a 10-fold cross-validation experiment on r1+∪r i−. in each experiment, we first used aleph, the ilp tool previously described, as a structure learner on the training set. with alchemy, 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 shows the averaged auc-roc and auc-pr values obtained within a large range of values of the hyperparameter λ of the mln while table 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. the table reports the averaged auc values and standard deviations obtained with mlns for thirty ten folds cross-validation experiments conducted on a regulatory network between the genes in ga. the results are reported for different values of the regularization parameter λ. the table reports the averaged auc values and standard deviations obtained with svms for thirty ten folds cross-validation experiments conducted on a regulatory network between the genes in ga. the column “sum” shows the results when the pairwise kernel is derived from the sum of genomic kernels, while the column “pairwise sum” shows the results obtained using the sum of pairwise kernels derived from each genomic kernel. the results are reported for different values of the regularization parameter c. 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 r2+ of regulations between the genes of interest on this date. we noticed that 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 we used the dataset r1+ from containing regulations as positive training set and tried to infer the new regulations in r2+∖r1+ using asymmetric bagging. to that end, we randomly sampled negative examples training sets r i−, i =  …, with r i−⊆r1+¯∖r2+ and |r i−|=|r1+|. 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/ of r1+ and r i− were considered for the training set and the remaining regulations were considered for the validation set. we computed each time the f1-measure obtained on the validation set for different threshold values between and 1: () f1= precision∗recallprecision+recall. we selected the threshold maximizing the averaged f1-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 new regulations. we averaged the predictions obtained and used the selected threshold value to compute the true positive rate. the tpr values obtained using bagged mln and bagged pairwise svms respectively are given in table 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. this table lists the true positive rates obtained for the prediction of regulations in r2+∖r1+ from r using asymmetric bagging with samples for bagged mlns and bagged pairwise svms. the tpr values were obtained using a threshold maximizing the averaged f1-measure on a validation set. notations are given in tables and 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. high confidence differentially expressed genes in prcid versus the corresponding control cells were identified. from these genes, we selected genes that were not part of ga 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 ga or with one gene of this new set. from these genes, we selected a subset of genes, called gb, that had at least a biological process annotation from go in common with genes from ga. 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 ga and gb, 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 r3+. we used the set r2+ of regulations from the dataset ga as the positive examples training set and used asymmetric bagging with samples on r for each sampling we applied the predictor on the sets r3+ and r3+¯, using descriptors on both set of genes ga and gb. we obtained score predictions for each interaction between one gene of ga and one gene of gb. table 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 this table lists the auc-roc and auc-pr values obtained for the prediction of regulations betweenga and gb for bagged mlns and bagged pairwise svms, with notations given in tables and these results were obtained using asymmetric bagging with samples on the set r 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 and negative examples 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. 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 and the implementation provided by the r package proc. the obtained p-values are given in the additional file we observe from this results that the p-values are less than 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 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 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 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. this table presents the results obtained with bagged mlns and bagged pairwise svms on the third task when using only gene expression data as gene descriptors. to conclude, we have shown in this section that bagged sum svm performs well in task and task while bagged pairwise sum svm performs well in task 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 and 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 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: protloccell ∧ expsiid2 ∧ expsiid2 ⇒ regulates processbio ∧ processbio ⇒ regulates expsiid2 ∧ expprcid2 ∧ expsiid2 ∧ expprcid2 ⇒ regulates expprcid2 ∧ expwt ∧ expprcid2 ⇒ regulates the first rule means that a gene overexpressed in transient knock-down of id 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 id compared to transient knock-down of id indicates a regulation between g and g regarding the last rule, it indicates that a high expression value of g in the prcid condition and the increase of the expression of g between wild-type condition and prcid 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.  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. on top of that, multiple kernel learning allows to select sources of informations thought the learning of sparse linear combination of kernels. 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, or using an extension of first-order logic to sequence variables like those of. 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 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. 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. finally the biologist interested in the id 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. competing interests the authors declare that they have no competing interests. authors’ contributions cb checked the encoding of the data in first-order logic, built the new datasets, developed codes for the whole pipeline and tested all the methods. fab and cv conceived the methodological part of the project, mad conceived its biological part about the id network. dc and mad did the experiments in the wet laboratory and provide a preliminary analysis of the data. jd encoded the data and the background knowledge into first-order logic. cv supervised the use of ilp for rule learning and fab supervised the mln approach, the asymmetric bagging part and designed with cb the baseline svm. mad participated to design the performance assessment with cb and fab. cb and fab drafted the document with the help of cv and mad. all authors read and approved the final manuscript. supplementary material additional file 1 list of gene symbols used in the study. click here for file additional file 2 choice of aleph parameters. click here for file additional file 3 roc and pr curves obtained for the prediction of regulations between setsga and gb. the figure represents the roc and pr curves obtained with bagged mlns, bagged svms and bagged svms using the hyperparameters associated with the best auc-roc values. click here for file additional file 4 p-values of the non-parametric test based on mann whitney statistics to compare auc-roc obtained with bagged mlns and bagged pairwise svms. click here for file acknowledgements we thank xavier gidrol and vincent frouin for their fruitful comments at the beginning of the study. this work was supported by the agence nationale de la recherche. the work of cb and fab was completed using funding from.
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. 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. once learnt, the pairwise classifier can be used to predict new regulations. in this work, we explore the framework of markov logic networks that combine features of probabilistic graphical models with the expressivity of first-order logic rules. 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. 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. 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. we propose to learn a markov logic network, e.g. finally our main result concerns a third test that measures the ability of the method to predict regulations with a new set of genes.
 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 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 that combine features of probabilistic graphical models with the expressivity of first-order logic rules.  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.  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. 