Abstract
Background
The regulation of gene expression by transcription factors is a key determinant of cellular phenotypes. Deciphering genomewide networks that capture which transcription factors regulate which genes is one of the major efforts towards understanding and accurate modeling of living systems. However, reverseengineering the network from gene expression profiles remains a challenge, because the data are noisy, high dimensional and sparse, and the regulation is often obscured by indirect connections.
Results
We introduce a gene regulatory network inference algorithm ENNET, which reverseengineers networks of transcriptional regulation from a variety of expression profiles with a superior accuracy compared to the stateoftheart methods. The proposed method relies on the boosting of regression stumps combined with a relative variable importance measure for the initial scoring of transcription factors with respect to each gene. Then, we propose a technique for using a distribution of the initial scores and information about knockouts to refine the predictions. We evaluated the proposed method on the DREAM3, DREAM4 and DREAM5 data sets and achieved higher accuracy than the winners of those competitions and other established methods.
Conclusions
Superior accuracy achieved on the three different benchmark data sets shows that ENNET is a top contender in the task of network inference. It is a versatile method that uses information about which gene was knockedout in which experiment if it is available, but remains the top performer even without such information. ENNET is available for download from https://github.com/slawekj/ennet webcite under the GNU GPLv3 license.
Keywords:
Gene regulatory networks; Network inference; Ensemble learning; BoostingBackground
Regulation of gene expression is a key driver of adaptation of living systems to changes in the environment and to external stimuli. Abnormalities in this highly coordinated process underlie many pathologies. At the transcription level, the control of the amount of mRNA transcripts involves epigenetic factors such as DNA methylation and, in eukaryotes, chromatin remodeling. But the key role in both prokaryotes and eukaryotes is played by transcription factors (TF), that is, proteins that can bind to DNA in the regulatory regions of specific genes and act as repressors or inducers of their expression. Many interactions between transcription factors and genes they regulate have been discovered through traditional molecular biology experiments. With the introduction of highthroughput experimental techniques for measuring gene expression, such as DNA microarrays and RNASeq, the goal moved to reverseengineering genomewide gene regulatory networks (GRNs) [1]. Knowledge of GRNs can facilitate finding mechanistic hypotheses about differences between phenotypes and sources of pathologies, and can help in the drug discovery and bioengineering.
High throughput techniques allow for collecting genomewide snapshots of gene expression across different experiments, such as diverse treatments or other perturbations to cells [2]. Analyzing these data to infer the regulatory network is one of the key challenges in the computational systems biology. The difficulty of this task arises from the nature of the data: they are typically noisy, high dimensional, and sparse [3]. Moreover, discovering direct causal relationships between genes in the presence of multiple indirect ones is not a trivial task given the limited number of knockouts and other controlled experiments. Attempts to solve this problem are motivated from a variety of different perspectives. Most existing computational methods are examples of influence modeling, where the expression of a target transcript is modeled as a function of the expression levels of some selected transcription factors. Such a model does not aim to describe physical interactions between molecules, but instead uses inductive reasoning to find a network of dependencies that could explain the regularities observed among the expression data. In other words, it does not explain mechanistically how transcription factors interact with regulated genes, but indicate candidate interactions with a strong evidence in expression data. This knowledge is crucial to prioritize detailed studies of the mechanics of the transcriptional regulation.
One group of existing methods describes GRN as a system of ordinary differential equations. The rate of change in expression of a transcript is given by a function of the concentration levels of transcription factors that regulate it. Network inference includes two steps: a selection of a model and an estimation of its parameters. Popular models imply linear functions a priori [47]. Bayesian Best Subset Regression (BBSR) [8] has been proposed as a novel model selection approach, which uses Bayesian Information Criterion (BIC) to select an optimal model for each target gene. Another group of methods employ probabilistic graphical models that analyze multivariate joint probability distributions over the observations, usually with the use of Bayesian Networks (BN) [911], or Markov Networks (MN) [12]. Various heuristic search schemes have been proposed in order to find parameters of the model, such as greedyhill climbing or the Markov Chain Monte Carlo approach [13]. However, because learning optimal Bayesian networks from expression data is computationally intensive, it remains impractical for genomewide networks.
Other approaches are motivated from statistics and information theory. TwixTwir [14] uses double twoway ttest to score transcriptional regulations. The nullmutant zscore algorithm [15] scores interactions based on a zscore transformed knockout expression matrix. Various algorithms rely on estimating and analyzing crosscorrelation and mutual information (MI) of gene expression in order to construct a GRN [1620], including ANOVA η^{2} method [21]. Improvements aimed at removing indirect edges from triples of genes have been proposed, including techniques such as the Data Processing Inequality in ARACNE [22,23], and the adaptive background correction in CLR [24]. Another method, NARROMI [25], eliminates redundant interactions from the MI matrix by applying ODEbased recursive optimization, which involves solving a standard linear programming model.
Recently, machinelearning theory has been used to formulate the network inference problem as a series of supervised gene selection procedures, where each gene in turn is designated as the target output. One example is MRNET [26], which applies the maximum relevance/minimum redundancy (MRMR) [27] principle to rank the set of transcription factors according to the difference between mutual information with the target transcript (maximum relevance) and the average mutual information with all the previously ranked transcription factors (minimum redundancy). GENIE3 [28] employs Random Forest algorithm to score important transcription factors, utilizing the embedded relative importance measure of input variables as a ranking criterion. TIGRESS [29] follows a similar approach but is based on the least angle regression (LARS). Recently, boosting [30,31] was also used to score the importance of transcription factors, in ADANET [32] and OKVARBoost [33] methods.
In this paper, we propose a method that combines gradient boosting with regression stumps, augmented with statistical reestimation procedures for prioritizing a selected subset of edges based on results from the machinelearning models. We evaluated our method on the DREAM3, DREAM4 and DREAM5 network inference data sets, and achieved results that in all cases were better than the currently available methods.
Methods
The ENNET algorithm
Formulating the gene network inference problem
The proposed algorithm returns a directed graph of regulatory interactions between P genes in form of a weighted adjacency matrix , where v_{i,j} represents regulation of gene j by gene i. As an input, it takes gene expression data from a set of experiments, together with the metadata describing the conditions of the experiments, including which genes were knocked out. Usually, the raw expression data need to be preprocessed before any inference method could be applied to reverseengineer a GRN. Preprocessing has a range of meanings, here it is regarded as a process of reducing variations or artifacts, which are not of the biological origin. It is especially important when the expression is measured with multiple highdensity microarrays [34]. Concentration levels of transcripts must be adjusted and the entire distribution of adjusted values aligned with a normal distribution. Methods for normalization of expression data are outside of the scope of our work. The data we used were already normalized using RMA [34,35] by the DREAM challenge organizers. We further normalized the expression data to zero mean and unit standard deviation.
The network inference process relies heavily on the type of expression data provided as an input. Two main groups of expression profiles are: the one with known, and the one with unknown initial perturbation state of the expression of genes in the underlying network of regulatory interactions. For example, knockout and knockdown data are provided with the additional metadata, which describe which genes were initially perturbed in each experiment. On the other hand, multifactorial and time series data are usually expression profiles of an unknown initial state of genes. Wildtype, knockout, knockdown, and multifactorial data describe the expression of initially perturbed genes, which are however in a steady state at the time of measurement, whereas time series data describe the dynamics of the expression levels of initially perturbed genes. The types of data available in popular benchmark data sets are summarized in Table 1.
Table 1. Different types of expression data provided in popular data sets
The variability of possible input scenarios poses a problem of representing and analyzing expression data. Here, we operate on an N×P expression matrix E, where e_{i,j} is the expression value of the jth gene in the ith sample. Columns of matrix E correspond to genes, rows correspond to experiments. We also define a binary perturbation matrix K, where k_{i,j} is a binary value corresponding to the jth gene in the ith sample, just like in the matrix E. If k_{i,j} is equal to 1, it means that the jth gene is known to be initially perturbed, for example knocked out, in the ith experiment. Otherwise k_{i,j} is equal to 0. If no information is available about knockouts, all values are set to 0.
Decomposing the inference problem into gene selection problems
We decompose the problem of inferring the network of regulatory interactions targeting all P genes into P independent subproblems. In each subproblem incoming edges from transcription factors to a single gene transcript are discovered. For the kth decomposed subproblem we create a target expression vector Y_{k} and a feature expression matrix X_{−k}. Columns of the X_{−k} matrix constitute a set of possible transcription factors. Vector Y_{k} corresponds to the expression of the transcript, which is possibly regulated by transcription factors from X_{−k}. In a single gene selection problem we decide which TFs contribute to the target gene expression across all the valid experiments. Columns of X_{−k} correspond to all the possible TFs, but if a target gene k is also a transcription factor, it is excluded from X_{−k}. We do not consider a situation in which a transcription factor would have a regulatory interaction with itself. When building the target vector Y_{k} corresponding to the kth target gene, k∈{1,...,P}, we consider all the experiments valid except from the ones in which the kth gene was initially perturbed, as specified in the perturbation matrix K. We reason that the expression value of the kth gene in those experiments is not determined by its TFs, but by the external perturbation. Each row in the Y_{k} vector is aligned with a corresponding row in the X_{−k} matrix. In order to justify all the possible interactions we need to solve a gene selection problem for each target gene. For example, if a regulatory network consists of four genes (P=4), we need to solve four gene selection problems. In the kth problem, k∈{1,2,3,4}, we find which TFs regulate the kth target gene. In other words, we calculate the kth column of the output adjacency matrix V.
Solving the gene selection problems
Once the target gene expression vector Y_{k} and the TF expression matrix X_{−k} are created for each gene k, we solve each kth gene selection problem independently, in the following way. We search for the subset of columns in X_{−k} that are related to the target vector Y_{k} by an unknown function f_{k}, as shown in Equation 1,
where ε_{k} is a random noise. A function f_{k} represents a pattern of regulatory interactions that drive the expression of the kth gene. We want f_{k} to rely only on a small number of genes acting as transcription factors, those that are the true regulators of gene k. Essentially, this is a feature selection or a gene selection task [28,32,36,37], where the goal is to model the target response Y_{k} with an optimal small set of important predictor variables, i.e., a subset of columns of the X_{−k} matrix. A more relaxed objective of the gene selection is the variable ranking, where the relative relevance for all input columns of the X_{−k} matrix is obtained with respect to the target vector Y_{k}. The higher a specific column is in that ranking, the higher the confidence that a corresponding TF is in a regulatory interaction with the target gene k.
Our solution to the variable ranking involves ensemble learning. We use an iterative regression method, which in each iteration chooses one transcription factor based on an optimality criterion, and adds it to the nonlinear regression ensemble. The main body of our method, presented in Figure 1, is based on Gradient Boosting Machine [38] with a squared error loss function. First, ENNET initializes f_{0} to be an optimal constant model, without selecting any transcription factor. In other words, f_{0} is initialized to an average of Y_{k}. At each next tth step the algorithm creates an updated model f_{t}, by fitting a base learner h_{t} and adding it to the previous model f_{t−1}. The base learner is fitted to a sample of pseudo residuals, with respect to a sample of transcription factors, and thus is expected to reduce the error of the model. Pseudoresiduals are recalculated at the beginning of each iteration with respect to the current approximation f_{t}. As a base learner, we use regression stumps, which select a single TF that best fits pseudo residuals. A regression stump h_{t}(x) partitions the expression values x of a candidate TF into two disjoint regions R_{1t} and R_{2t}, where , and returns values γ_{1t} and γ_{2t}, respectively, for those regions, as shown in Equation 2,
Figure 1. The flowchart of the ENNET algorithm. ENNET algorithm is a modification of a Gradient Boosting Machine algorithm, with a squared error loss function and a regression stump base learner. The algorithm calculates a vector of importance scores of transcription factors, which can possibly regulate a target gene. It is invoked P times in a problem of inferring a Pgene network, i.e., a Pcolumn adjacency matrix V.
where I is the identity function returning the numerical 1 for the logical true, and the numerical 0 for the logical false. Regions R_{1t}, R_{2t} are induced such that the leastsquares improvement criterion is maximized:
where w_{1t}, w_{2t} are proportional to the number of observations in regions R_{1t}, R_{2t} respectively, and γ_{1t}, γ_{2t} are corresponding response means. That is, γ_{1t} is the average of the values from the vector of pseudoresiduals for those samples where an expression of the chosen TF falls into the region R_{1t}. The value of γ_{2t} is defined in an analogous way. The averages γ_{1t} and γ_{2t} are used as the regression output values for regions R_{1t} and R_{2t}, respectively, as shown in Equation 2. The criterion in Equation 3 is evaluated for each TF, and the transcription factor with the highest improvement is selected. In each tth step, we only use a random portion of rows and columns of X_{−k}, sampled according to the observation sampling rate s_{s}, and the TF sampling rate s_{f}.
The procedure outlined above creates a nonlinear regression model of the target gene expression based on the expression of transcription factors. However, in the network inference, we are interested not in the regression model as a whole, but only in the selected transcription factors. In each tth step of the ENNET algorithm, only one TF is selected as the optimal predictor. The details of the regression model can be used to rank the selected TFs by their importance. Specifically, if a transcription factor φ_{t} is selected in an iteration t, an improvement serves as an importance score for that φ_{t}th TF. If the same TF is selected multiple times at different iterations, its final importance score is a sum of the individual scores.
In the training of the regression model, the parameter ν, known as a shrinkage factor, is used to scale a contribution of each tree by a factor ν∈(0,1) when it is added to the current approximation. In other words, ν controls the learning rate of the boosting procedure. Shrinkage techniques are also commonly used in neural networks. Smaller values of ν result in a larger training risk for the same number of iterations T. However, it has been found [38] that smaller values of ν reduce the test error, and require correspondingly larger values of T, which results in a higher computational overhead. There is a tradeoff between these two parameters.
Refining the inferred network
Once the solutions of the independent gene selection problems are calculated, we compose the adjacency matrix V representing a graph of inferred regulatory interactions. Each of the solutions constitutes a single columnvector, therefore we obtain the adjacency matrix V by binding all the partial solutions columnwise. Then we apply a reevaluation algorithm to achieve an improved final result. The first step does not require any additional data to operate other than the previously calculated adjacency matrix V. It exploits the variance of edge probabilities in the rows of V, i.e., edges outgoing from a single transcription factor, as a measure of the effect of transcriptional regulation. We score transcription factors based on their effects on multiple targets. We assume that the effect of transcriptional regulation on a directly regulated transcript is stronger than the one of the regulation on indirectly regulated transcripts, e.g. transcripts regulated through another transcription factor. Otherwise, knocking out a single gene in a strongly connected component in a network of regulatory interactions would cause the same rate of perturbation of the expression level of all the transcripts in that component. As a measure of that effect we use previously a calculated adjacency matrix V and multiply each row of V matrix by its variance . An updated adjacency matrix V^{1} is given by Equation 4:
where is a variance in the ith row of V. Note that V matrix is built columnwise, i.e., a single column of V contains the relative importance scores of all the transcription factors averaged over all the base learners with respect to a single target transcript. On the other hand, rows of V matrix are calculated independently in different subproblems of the proposed inference method. Each row of V contains relative importance scores with respect to a different target transcript. We reason that if a transcription factor regulates many target transcripts, e.g. a transcription factor is a hub node, the variance in a row of V corresponding to that transcription factor is elevated and therefore it indicates an important transcription factor.
The second step of refining the network requires knockout expression data. We reason that direct regulation of a transcript by a transcription factor would lead to a distinct signature in the expression data if the transcription factor was knocked out. A similar reasoning gave foundations for the nullmutant zscore method [15] of reverseengineering GRNs. However, in the proposed method this step is only applied if knockout expression profiles are available. In this step we calculate an adjacency matrix V^{2}, which is an update to an already derived adjacency matrix V^{1}, as shown in Equation 5:
where is an average expression value of the jth transcript in all the experiments α(i) in which the ith gene was knockedout, as defined by K matrix, is the mean expression value for that transcript across all the other knockout experiments, β(i), and σ_{j} is the standard deviation of the expression value of that transcript in all the knockout experiments. The coefficient shows how many standard deviations the typical expression of the jth transcript was different from the average expression in the experiment in which its potential ith transcription factor was knockedout.
Performance evaluation
A considerable attention has been devoted in recent years to the problem of evaluating performance of the inference methods on adequate benchmarks [35,39]. The most popular benchmarks are derived from wellstudied in vivo networks of model organisms, such as E. coli[40] and S. cerevisiae[41], as well as artificially simulated in silico networks [39,4245]. The main disadvantage of in vivo benchmark networks is the fact that experimentally confirmed pathways can never be assumed complete, regardless of how well the model organism is studied. Such networks are assembled from known transcriptional interactions with strong experimental support. As a consequence, gold standard networks are expected to have few false positives. However, they contain only a subset of the true interactions, i.e., they are likely to contain many false negatives. For this reason, artificially simulated in silico networks are most commonly used to evaluate network inference methods. Simulators [39] mimic real biological systems in terms of topological properties observed in biological in vivo networks, such as modularity [46] and occurrences of network motifs [47]. They are also endowed with dynamical models of a transcriptional regulation, thanks to the use of nonlinear differential equations and other approaches [42,48,49], and consider both transcription and translation processes in their dynamical models [4850] using a thermodynamic approach. Expression data can be generated deterministically or stochastically and experimental noise, such as the one observed in microarrays, can be added [51].
Here, we used several popular benchmark GRNs to evaluate the accuracy of our proposed algorithm and compare it with the other inference methods. The data sets we used come from Dialogue for Reverse Engineering Assessments and Methods (DREAM) challenges and are summarized in Table 1. We evaluated the accuracy of the methods using the Overall Score metric proposed by the authors of DREAM challenges [35], as shown in Equation 6:
where and are geometric means of pvalues of networks constituting each DREAM challenge, relating to an area under the PrecisionRecall curve (AUPR) and an area under the ROC curve (AUROC), respectively.
Results and discussion
We assessed the performance of the proposed inference algorithm on large, universally recognized benchmark networks of 100 and more genes, and compared it to the stateoftheart methods. We summarize the results of running different inference methods in Figure 2. For a comparison we selected a range of established methods from literature: ARACNE, CLR, and MRNET as implemented in the minet R package [52], GENIE3 and C3NET as implemented by their respective authors, our previously reported method ADANET, and the top three performers in each of the three DREAM challenges as listed on the DREAM web site. Some of the methods were designed for use with knockout data, while others are developed with multifactorial data in mind, where no information is given about the nature of the perturbations. Therefore, depending on the nature of the particular DREAM data set, only the suitable group of methods is used for the comparison.
Figure 2. The Overall Score of GRN inference methods by data set. Results of the different inference methods on DREAM challenges. Results of the stateoftheart methods were collected after running the algorithms with the default sets of parameters on preprocessed data. Results in the “Winner of the challenge” part of the figure correspond to the best methods participating in the challenge.
The accuracy of ENNET
DREAM3 [15,53,54] features in silico networks and expression data simulated using GeneNetWeaver software. Benchmark networks were derived as subnetworks of a system of regulatory interactions from known model organisms: E. coli and S. cerevisiae. In this study we focus on a DREAM3 size 100 subchallenge, as the largest of DREAM3 suite. The results of all the competing methods except those that are aimed at multifactorial problems are summarized in Table 2. ENNET and Yip et al. methods achieved the best Overall Scores for that subchallenge, as well as the best scores for all the individual networks. However, it is believed from the analysis of the later challenges [39] that Yip et al. method made a strong assumption on the Gaussian type of a measurement noise, which was used in DREAM3, but was no longer used in later DREAM challenges. For example, in DREAM4 challenge Yip et al. method was ranked 7th.
Table 2. Results of the different inference methods on DREAM3 networks, challenge size 100
DREAM4 challenge [15,53,54] was posted one year after DREAM3 challenge. It features two large subchallenges: DREAM4 size 100, and DREAM4 size 100 multifactorial. For each subchallenge, the topology of the benchmark networks were derived from the transcriptional regulatory system of E. coli and S. cerevisiae. In DREAM4 size 100 subchallenge all the data types listed in Table 1 were available except multifactorial, therefore ADANET, GENIE3, CLR, C3NET, MRNET, and ARACNE methods were excluded from the comparison. The results of all the methods are summarized in Table 3. ENNET method clearly outperformed all the others and achieved consistently high scores across all the benchmark networks. In the second DREAM4 large subchallenge, DREAM4 size 100 multifactorial, only multifactorial data were available, therefore all the methods were included in the comparison, and run as originally designed. The results of all the methods are summarized in Table 4. ENNET achieved the best Overall Score.
Table 3. Results of the different inference methods on DREAM4 networks, challenge size 100
Table 4. Results of the different inference methods on DREAM4 networks, challenge size 100 multifactorial
Three benchmark networks in DREAM5 [35] were different in size, and structured with respect to different model organisms. However, this time expression data of the only one network were simulated in silico, the two other sets of expression data were measured in real experiments in vivo. Like in all DREAM challenges, in silico expression data were simulated using an opensource GeneNetWeaver simulator [54]. However, DREAM5 was the first challenge where participants were asked to infer GRNs on a genomic scale, e.g. for thousands of target genes, and hundreds of known transcription factors. Gold standard networks were obtained from two sources: RegulonDB database [40], and Gene Ontology (GO) annotations [55]. The results of all the inference methods for DREAM5 expression data are summarized in Table 5. ENNET achieved the best score for the in silico network, and the best Overall Score, as well as the best individual AUROC scores for all the networks. Clearly all the participating methods achieved better scores for an in silico network than for either one of in vivo networks. ENNET shows better in vivo results than the other methods in terms of an area under the the ROC curve. Still, predictions for in vivo expression profiles show a low overall accuracy. One of the reasons for a poor performance of the inference methods for such expression profiles is a fact that experimentally confirmed pathways, and consequently gold standards derived from them, cannot be assumed complete, regardless of how well is a model organism known. Additionally, there are regulators of gene expression other than transcription factors, such as miRNA, and siRNA. As shown in this study, in silico expression profiles provide enough information to confidently reverseengineer their underlying structure, whereas in vivo data hide a much more complex system of regulatory interactions.
Table 5. Results of the different inference methods on DREAM5 networks
Computational complexity of ENNET
Computational complexity of ENNET depends mainly on the computational complexity of the regression stump base learner, which is used in the main loop of the algorithm. As shown in Figure 1, we call the regression stump algorithm T times for each kth target gene, k∈{1,...,P}. Given a sorted input, a regression stump is O(PN) complex. We sort the expression matrix in an O(PN logN) time. All the other instructions in the main loop of ENNET are at most O(N). The computational complexity of the whole method is thus O(PN logN+TP^{2}N+TPN). Because, in practice, the dominating part of the sum is TP^{2}N, we report a final computational complexity of ENNET as O(TP^{2}N), and compare it to the other inference methods in Table 6. Note that the measure for the informationtheoretic methods: CLR, MRNET, and ARACNE does not include a calculation of the mutual information matrix.
Table 6. The computational complexity of ENNET and the other GRN inference methods
When implementing ENNET algorithm we took advantage of the fact that gene selection problems are independent of each other. Our implementation of the algorithm is able to calculate them in parallel if multiple processing units are available. User can choose from variety of parallel backends including multicore package for a single computer and parallelization based on Message Passing Interface for a cluster of computers. The biggest data we provided as input in our tests were in vivo expression profiles of S. cerevisiae from the DREAM 5 challenge. These are genomewide expression profiles of 5950 genes (333 of them are known transcription factors) measured in 536 experiments. It took 113 minutes and 30 seconds to calculate the network on a standard desktop workstation with one Intel®;Core™i7870 processor with 4 cores and two threads per core (in total 8 logical processors) and 16 GB RAM. However, it took only 16 minutes and 40 seconds to calculate the same network on a machine with four AMD Opteron™6282 SE processors, each with 8 cores and two threads per core (in total 64 logical processors) and 256 GB RAM. All the data sets from the DREAM 3 and the DREAM 4 challenges were considerably smaller, up to 100 genes. It took less than one minute to calculate each of these networks on a desktop machine.
Setting parameters of ENNET
The ENNET algorithm is controlled by four parameters: the two sampling rates s_{s} and s_{f}, the number of iterations T and the learning rate ν. The sampling rate of samples s_{s} and the sampling rate of transcription factors s_{f} govern the level of randomness when selecting, respectively, rows and columns of the expression matrix to fit a regression model. The default choice of the value of s_{s} is 1, i.e., we select with replacement a bootstrap sample of observations of the same size as an original training set at each iteration. Because some observations are selected more than once, around 0.37 of random training samples are out of bag in each iteration. It is more difficult to choose an optimal value of s_{f}, which governs how many transcription factors are used to fit each base learner. Setting this parameter to a low value forces ENNET to score transcription factors, even if their improvement criterion, as shown in Equation 2, would not have promoted them in a pure greedy search, i.e., s_{f}=1. However, if a chance of selecting a true transcription factor as a feature is too low, ENNET will suffer from selecting random genes as true regulators.
Even though reverseengineering of GRNs does not explicitly target a problem of predicting gene expression, we choose the values of sampling rates such that the squarederror loss of a prediction of the target gene expression as given by f_{T} (see Figure 1) is minimal. This is done without looking at the ground truth of regulatory connections. For each benchmark challenge we performed a grid search over (s_{s},s_{f})∈{0.1,0.3,0.5,0.7,1}×{0.1,0.3,0.5,0.7,1} with fixed ν=0.001, T=5000. For each specific set of parameters we analyzed an average 5fold crossvalidated loss over all the observations (across all gene selection problems). We further analyze our approach with respect to one of the challenges: DREAM4 size 100, as shown in Figure 3. The minimal average loss was achieved for s_{s}=1 and s_{f}=0.3 (see Figure 3 A), which is consistent with the default parameters proposed for Random Forest algorithm [28]. We also compared the measure based on an average loss with the Overall Score as defined by Equation 6. The results were consistent across the two measures, i.e., a selection of parameters that gave a low average loss also led to the accurate network predictions (see Figure 3 B). An advantage of the average loss measure is a fact that the gold standard network is not used to tune parameters.
Figure 3. The analysis of the sampling rates s_{s} and s_{f } for DREAM 4 size 100 challenge: the Overall Score and a loss. The analysis of the sampling rates s_{s} and s_{f} for the DREAM 4 size 100 challenge. A: For each set of parameters (s_{s}, s_{f}, M, ν) ∈ {0.1,0.3,0.5,0.7,1} × {0.1,0.3,0.5,0.7,1} × {5000} × {0.001} we analyzed an average 5fold crossvalidated loss over all the observations (across all gene selection problems) from all 5 networks. The minimal average loss was achieved for high values of s_{s} = 1 and low values of s_{f} = 0.3. B: We also compared the measure based on an average loss with the original Overall Score, as proposed by the authors of the DREAM challenge. The results were consistent across the two measures, i.e., the parameters that gave low average loss also led to accurate network predictions (a high Overall Score).
In Figure 4 we present a detailed analysis of the accuracy of the GRN inference across different networks of the DREAM4 size 100 challenge. Each point on both Figure 4 A and Figure 4 B is a result of running ENNET with different parameters: (s_{s},s_{f})∈{0.1,0.3,0.5,0.7,1}×{0.1,0.3,0.5,0.7,1} with fixed ν=0.001, T=5000. The highlighted points are corresponding to s_{s}=1, s_{f}=0.3, ν=0.001, T=5000. An area under the PrecisionRecall curve and an area under the ROC curve are two different measures of the accuracy of an inferred network, which are well preserved across the five networks: for each separate network we observe that AUPR and AUROC decreases in a function of an average loss. As the Overall Score is closely related to AUPR and AUROC, the results shown in Figure 4 explain the shape of a surface shown in Figure 3.
Figure 4. The analysis of the sampling rates s_{s }and s_{f} for DREAM 4 size 100 challenge: AUPR, AUROC, and a loss. The analysis of the sampling rates s_{s} and s_{f} for DREAM 4 size 100 challenge. A: For each set of parameters (s_{s}, s_{f}, M, ν) ∈ {0.1,0.3,0.5,0.7,1} × {0.1,0.3,0.5,0.7,1} × {5000} × {0.001} we analyzed an area under the PrecisionRecall curve (AUPR) in function of an average 5fold crossvalidated loss over all the observations (across all gene selection problems) from all 5 networks. For each network AUPR is decreasing in a function of a loss. For each network a point corresponding to the default set of parameters is highlighted, i.e., (s_{s}, s_{f}, M, ν) = (1,0.3,5000,0.001). Usually, the default set of parameters gives the minimal loss (maximal AUPR). B: By analogy, different choices of parameters lead to a different area under the ROC curve (AUROC). The two measures are consistent with each other.
As ENNET uses boosting, it needs a careful tuning of the number of iterations T and the learning rate ν. It has been shown [38] that parameters T and ν are closely coupled. Usually the best prediction results are achieved when ν is fixed to a small positive number, e.g. ν≤0.001, and the optimal value of TY is found in a process of crossvalidation. As described above, we reason that the choice of parameters, which gives a low average loss on a crossvalidated test set, leads to an accurate network prediction. Therefore in Figure 5 we present how an average loss depends on T∈{1,...,5000} for different values of ν∈{0.001,0.005,0.01,0.05,0.1}, with fixed s_{s}=1, s_{f}=0.3. Each of the line shows how much ENNET overtrains the data for a given T and ν. Finally, the optimal choice of parameters for DREAM4 size 100 challenge is s_{s}=1, s_{f}=0.3, T=5000, ν=0.001. Following the same practice, we used this default set of parameters: s_{s}=1, s_{f}=0.3, T=5000, ν=0.001 to evaluate ENNET algorithm on all the benchmark networks using ground truth, i.e., for calculating the Overall Score and comparing it to the other algorithms.
Figure 5. The analysis of the number of iterations T and the shrinkage ν for the DREAM 4 size 100 challenge. The analysis of the number of iterations T and the shrinkage factor ν for DREAM 4 size 100 challenge. These two parameters are closely coupled: the lower is the shrinkage parameter ν, the more iterations T are needed to train the model such that it achieves the minimal loss.
Stability of ENNET
Because ENNET uses random sampling of samples and features at each iteration of the main loop, as shown in Figure 1, it may calculate two different networks for two different executions on the same expression data. With the default choice of parameters, i.e., s_{s}=1, s_{f}=0.3, T=5000, ν=0.001, we expect numerous random resamplings, and therefore we need to know if a GRN calculated by ENNET is stable between different executions. We applied ENNET to the 5 networks that form DREAM 4 size 100 benchmark, repeating the inference calculations independently ten times for each network. Then, for each network, we calculated a Spearman’s rank correlation between all pairs among the ten independent runs. The lowest correlation coefficient we obtained was ρ>0.975, with pvalue <2.2e−16, indicating that the networks that result from independent runs are very similar. This proves that ENNET, despite being a randomized algorithm, finds a stable solution to the inference problem.
Conclusions
We have proposed the ENNET algorithm for reverseengineering of Gene Regulatory Networks. ENNET uses a variety of types of expression data as an input, and shows robust performance across different benchmark networks. Moreover, it does not assume any specific model of a regulatory interaction and do not require finetuning of its parameters, i.e., we define the default set of parameters, which promises accurate predictions for the future networks. Nevertheless, together with the algorithm, we propose a procedure of tuning parameters of ENNET towards minimizing empirical loss. Processing genomescale expression profiles is feasible with ENNET: including up to a few hundred transcription factors, and up to a few thousand regulated genes. As shown in this study, the proposed method compares favorably to the stateoftheart algorithms on the universally recognized benchmark data sets.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
JS and TA conceived the method and drafted the manuscript. JS implemented the method and ran the experiments. JS and TA read and approved the final manuscript.
References

Someren E, Wessels L, Backer E, Reinders M: Genetic network modeling.
Pharmacogenomics 2002, 3(4):507525. PubMed Abstract  Publisher Full Text

Eisen M, Spellman P, Brown P, Botstein D: Cluster analysis and display of genomewide expression patterns.
Proc Natl Acad Sci 1998, 95(25):14863. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gardner T, Faith J: Reverseengineering transcription control networks.
Phys Life Rev 2005, 2:6588. PubMed Abstract  Publisher Full Text

Chen T, He H, Church G, et al.: Modeling gene expression with differential equations. In Pacific Symposium on Biocomputing, Volume 4. Singapore: World Scientific Press; 1999:44.

D’haeseleer P, Wen X, Fuhrman S, Somogyi R, et al.: Linear modeling of mRNA expression levels during CNS development and injury,. In Pacific Symposium on Biocomputing, Volume 4. Singapore: World Scientific Press; 1999:4152.

Gardner TS, di Bernardo D, Lorenz D, Collins JJ: Inferring genetic networks and identifying compound mode of action via expression profiling.
Science 2003, 301(5629):102105. PubMed Abstract  Publisher Full Text

Yip K, Alexander R, Yan K, Gerstein M: Improved reconstruction of in silico gene regulatory networks by integrating knockout and perturbation data.
PLoS One 2010, 5:e8121. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Greenfield A, Hafemeister C, Bonneau R: Robust datadriven incorporation of prior knowledge into the inference of dynamic regulatory networks.
Bioinformatics 2013, 29(8):10601067. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Friedman N, Linial M, Nachman I, Pe’er D: Using Bayesian networks to analyze expression data.
J Comput Biol 2000, 7(3–4):601620. PubMed Abstract  Publisher Full Text

Perrin BE, Ralaivola L, Mazurie A, Bottani S, Mallet J, d‘Alche Buc F: Gene networks inference using dynamic Bayesian networks.
Bioinformatics 2003, 19(suppl 2):ii138ii148. PubMed Abstract  Publisher Full Text

Yu J, Smith V, Wang P, Hartemink A, Jarvis E: Advances to Bayesian, network inference for generating causal networks from observational biological data.
Bioinformatics 2004, 20(18):35943603. PubMed Abstract  Publisher Full Text

Segal E, Wang H, Koller D: Discovering molecular pathways from protein interaction and gene expression data.
Bioinformatics 2003, 19(suppl 1):i264i272. PubMed Abstract  Publisher Full Text

Neapolitan R: Learning Bayesian Networks. Upper Saddle River: Pearson Prentice Hall; 2004.

Qi J, Michoel T: Contextspecific transcriptional regulatory network inference from global gene expression maps using double twoway ttests.
Bioinformatics 2012, 28(18):23252332. PubMed Abstract  Publisher Full Text

Prill R, Marbach D, SaezRodriguez J, Sorger P, Alexopoulos L, Xue X, Clarke N, AltanBonnet G, Stolovitzky G: Towards a rigorous assessment of systems biology models: the DREAM3 challenges.
PLoS One 2010, 5(2):e9202. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bansal M, Belcastro V, AmbesiImpiombato A, di Bernardo D: How to infer gene networks from expression profiles.
Mol Syst Biol 2007, 3:78. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Butte A, Kohane I: Mutual information relevance networks: functional genomic clustering using pairwise entropy measurements. In Pacific Symposium on Biocomputing, Volume 5. Singapore: World Scientific Press; 2000:418429.

Lee W, Tzou W: Computational methods for discovering gene networks from expression data.
Brief Bioinform 2009, 10(4):408423. PubMed Abstract  Publisher Full Text

Markowetz F, Spang R: Inferring cellular networks–a review.
BMC Bioinformatics 2007, 8(Suppl 6):S5. BioMed Central Full Text

Altay G, EmmertStreib F: Inferring the conservative causal core of gene regulatory networks.
BMC Syst Biol 2010, 4:132. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Küffner R, Petri T, Tavakkolkhah P, Windhager L, Zimmer R: Inferring gene regulatory networks by ANOVA.
Bioinformatics 2012, 28(10):13761382. PubMed Abstract  Publisher Full Text

Margolin A, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Favera R, Califano A: ARACNE: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context.
BMC Bioinformatics 2006, 7(Suppl 1):S7. BioMed Central Full Text

Margolin A, Wang K, Lim W, Kustagi M, Nemenman I, Califano A: Reverse engineering cellular networks.
Nat Protoc 2006, 1(2):662671. PubMed Abstract  Publisher Full Text

Faith J, Hayete B, Thaden J, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins J, Gardner T: Largescale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles.
PLoS Biol 2007, 5:e8. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zhang X, Liu K, Liu ZP, Duval B, Richer JM, Zhao XM, Hao JK, Chen L: NARROMI: a noise and redundancy reduction technique improves accuracy of gene regulatory network inference.
Bioinformatics 2013, 29:106113. PubMed Abstract  Publisher Full Text

Meyer P, Kontos K, Lafitte F, Bontempi G: Informationtheoretic inference of large transcriptional regulatory networks.

Ding C, Peng H: Minimum redundancy feature selection from microarray gene expression data. In Computational Systems Bioinformatics Conference CSB2003. Washington: IEEE; 2003:523528.

Irrthum A, Wehenkel L, Geurts P, et al.: Inferring regulatory networks from expression data using treebased methods.
PLoS One 2010, 5(9):e12776. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Haury AC, Mordelet F, VeraLicona P, Vert JP: TIGRESS: trustful inference of gene regulation using stability selection.
BMC Syst Biol 2012, 6:145. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Freund Y, Schapire RE: Experiments with a new boosting algorithm.

Freund Y, Schapire RE: A decisiontheoretic generalization of online learning and an application to boosting.
J Comput Syst Sci 1997, 55:119139. Publisher Full Text

Sławek J, Arodź T: ADANET: inferring gene regulatory networks using ensemble classifiers. In Proceedings of the ACM Conference on Bioinformatics, Computational Biology and Biomedicine. New York: ACM; 2012:434441.

Lim N, Şenbabaoğlu Y, Michailidis G, d’Alché Buc F: OKVARBoost: a novel boosting algorithm to infer nonlinear dynamics and interactions in gene regulatory networks.
Bioinformatics 2013, 29(11):14161423. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bolstad B, Irizarry R, Åstrand M Speed T: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias.
Bioinformatics 2003, 19(2):185193. PubMed Abstract  Publisher Full Text

Marbach D, Costello J, Küffner R, Vega N, Prill R, Camacho D, Allison K, Kellis M, Collins J, Stolovitzky G, et al.: Wisdom of crowds for robust gene network inference.

Theodoridis S, Koutroumbas K: Pattern Recognition. London: Elsevier/Academic Press; 2006.

Tuv E, Borisov A, Runger G, Torkkola K: Feature selection with ensembles, artificial variables, and redundancy elimination.

Friedman JH: Greedy function approximation: a gradient boosting machine.

Schaffter T, Marbach D, Floreano D: GeneNetWeaver: in silico benchmark generation and performance profiling of network inference methods.
Bioinformatics 2011, 27(16):22632270. PubMed Abstract  Publisher Full Text

GamaCastro S, Salgado H, PeraltaGil M, SantosZavaleta A, MuñizRascado L, SolanoLira H, JimenezJacinto V, Weiss V, GarcíaSotelo J, LópezFuentes A, et al.: RegulonDB version 7.0: transcriptional regulation of Escherichia coli K12 integrated within genetic sensory response units (gensor units).
Nucleic Acids Res 2011, 39(suppl 1):D98D105. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kim S, Imoto S, Miyano S: Inferring gene networks from time series microarray data using dynamic Bayesian networks.
Brief Bioinform 2003, 4(3):228235. PubMed Abstract  Publisher Full Text

Di Camillo B, Toffolo G, Cobelli C: A gene network simulator to assess reverse engineering algorithms.
Ann N Y Acad Sci 2009, 1158:125142. PubMed Abstract  Publisher Full Text

Kremling A Fischer S, Gadkar K, Doyle F, Sauter T, Bullinger E, Allgöwer F, Gilles E: A benchmark for fethods in reverse engineering and model discrimination: problem formulation and solutions.
Genome Res 2004, 14(9):17731785. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mendes P, Sha W, Ye K: Artificial gene networks for objective comparison of analysis algorithms.
Bioinformatics 2003, 19(suppl 2):ii122ii129. PubMed Abstract  Publisher Full Text

Van den Bulcke T, Van Leemput K, Naudts B, Van Remortel P, Ma H, Verschoren A, De Moor B, Marchal K: SynTReN a generator of synthetic gene expression data for design and analysis of structure learning algorithms.
BMC Bioinformatics 2006, 7:43. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Ravasz E, Somera A, Mongru D, Oltvai Z, Barabási A: Hierarchical organization of modularity in metabolic networks.
Science 2002, 297(5586):15511555. PubMed Abstract  Publisher Full Text

ShenOrr S, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli.
Nat Genet 2002, 31:6468. PubMed Abstract  Publisher Full Text

Hache H, Wierling C, Lehrach H, Herwig R: GeNGe: systematic generation of gene regulatory networks.
Bioinformatics 2009, 25(9):12051207. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Roy S, WernerWashburne M, Lane T: A system for generating transcription regulatory networks with combinatorial control of transcription.
Bioinformatics 2008, 24(10):13181320. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Haynes B, Brent M: Benchmarking regulatory network reconstruction with GRENDEL.
Bioinformatics 2009, 25(6):801807. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stolovitzky G, Kundaje A, Held G, Duggar K, Haudenschild C, Zhou D, Vasicek T, Smith K, Aderem A, Roach J: Statistical analysis of MPSS, measurements: application to the study of LPSactivated macrophage gene expression.
Proc Natl Acad Sci USA 2005, 102(5):14021407. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Meyer PE, Lafitte F, Bontempi G: Minet: AR/Bioconductor package for inferring large transcriptional networks using mutual information.
BMC Bioinformatics 2008, 9:461. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Marbach D, Prill R, Schaffter T, Mattiussi C, Floreano D, Stolovitzky G: Revealing strengths and weaknesses of methods for gene network inference.
Proc Natl Acad Sci 2010, 107(14):62866291. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Marbach D, Schaffter T, Mattiussi C, Floreano D: Generating realistic in silico gene networks for performance assessment of reverse engineering methods.
J Comput Biol 2009, 16(2):229239. PubMed Abstract  Publisher Full Text

Ashburner M, Ball C, Blake J, Botstein D, Butler H, Cherry J, Davis A, Dolinski K, Dwight S, Eppig J, et al.: Gene Ontology: tool for the unification of biology.
Nat Genet 2000, 25:25. PubMed Abstract  Publisher Full Text  PubMed Central Full Text