 Methodology article
 Open Access
 Published:
TIGRESS: Trustful Inference of Gene REgulation using Stability Selection
BMC Systems Biology volume 6, Article number: 145 (2012)
Abstract
Background
Inferring the structure of gene regulatory networks (GRN) from a collection of gene expression data has many potential applications, from the elucidation of complex biological processes to the identification of potential drug targets. It is however a notoriously difficult problem, for which the many existing methods reach limited accuracy.
Results
In this paper, we formulate GRN inference as a sparse regression problem and investigate the performance of a popular feature selection method, least angle regression (LARS) combined with stability selection, for that purpose. We introduce a novel, robust and accurate scoring technique for stability selection, which improves the performance of feature selection with LARS. The resulting method, which we call TIGRESS (for Trustful Inference of Gene REgulation with Stability Selection), was ranked among the top GRN inference methods in the DREAM5 gene network inference challenge. In particular, TIGRESS was evaluated to be the best linear regressionbased method in the challenge. We investigate in depth the influence of the various parameters of the method, and show that a fine parameter tuning can lead to significant improvements and stateoftheart performance for GRN inference, in both directed and undirected settings.
Conclusions
TIGRESS reaches stateoftheart performance on benchmark data, including both in silico and in vivo (E. coli and S. cerevisiae) networks. This study confirms the potential of feature selection techniques for GRN inference. Code and data are available on http://cbio.ensmp.fr/tigress. Moreover, TIGRESS can be run online through the GenePattern platform (GPDREAM, http://dream.broadinstitute.org).
Background
In order to meet their needs and adapt to changing environments, cells have developed various mechanisms to regulate the production of the thousands of proteins they can synthesize. Among them, the regulation of gene expression by transcription factors (TF) is preponderant: by binding to the promoter regions of their target genes (TG), TF can activate or inhibit their expression. Deciphering and understanding TFTG interactions has many potential farreaching applications in biology and medicine, ranging from the in silico modeling and simulation of the gene regulatory network (GRN) to the identification of new potential drug targets. However, while many TFTG interactions have been experimentally characterized in model organisms, the systematic experimental characterization of the full GRN remains a daunting task due to the large number of potential regulations.
The development of highthroughput methods, in particular DNA microarrays, to monitor gene expression on a genomewide scale has promoted new strategies to elucidate GRN. By systematically assessing how gene expression varies in different experimental conditions, one can try to reverse engineer the TFTG interactions responsible for the observed variations. Not surprisingly, many different approaches have been proposed in the last decade to solve this GRN reverse engineering problem from collections of gene expression data. When expression data are collected over time, for example, several methods have been proposed to construct dynamic models where TFTG interactions dictate how the expression level of a TG at a given time allows to predict the expression levels of its TG in subsequent times [1–11]. When expression data are not limited to time series, many methods attempt to capture statistical association between the expression levels of TG and candidate TF using correlation or informationtheoretic measures such as mutual information [12–14] or take explicit advantage of perturbations in the experiments such as gene knockdowns [15]. The difficulty to separate direct from indirect regulations has been addressed with the formalism of Bayesian networks [16–19], or by formulating the GRN inference problem as a feature selection problem [20]. Mutual informationbased ARACNE [13] was also designed to eliminate redundant edges. We refer to [21, 22] for detailed reviews and comparisons of existing methods.
Recent benchmarks and challenges have highlighted the good performance of methods which formalize the GRN inference problem as a regression and feature selection problem, namely, identifying a small set of TF whose expression levels are sufficient to predict the expression level of each TG of interest. The general idea that edges in a directed graph can be discovered node by node was addressed in, e.g.,[23]. Regarding the GRN inference application, this idea underlies the Bayesian network formalism [16], but is more directly addressed by GENIE3 [20], a method which uses random forests to identify TF whose expression levels are predictive for the expression level of each TG, and which is now recognized as stateoftheart on several benchmarks [20, 22]. Feature selection with random forests remains however poorly understood theoretically, and one may wonder how other wellestablished statistical and machine learning techniques for feature selection would behave to solve the GRN inference problem.
In this paper, we investigate the performance of a popular feature selection method, least angle regression (LARS) [24] combined with stability selection [25, 26], for GRN inference. LARS is a computationally efficient procedure for multivariate feature selection, closely related to Lasso regression [27]. Stability selection consists in running LARS or Lasso many times, resampling the samples and the variables at each run, and in computing the frequency with which each variable was selected across the runs. We introduce a novel, robust and accurate scoring technique for stability selection, which improves the performance of feature selection with LARS. The resulting method, which we call TIGRESS (for Trustful Inference of Gene REgulation with Stability Selection), was ranked among the top GRN inference methods in the DREAM5 gene reconstruction challenge and was evaluated to be the best linear regression based method [28]. We furthermore investigate in depth the influence of the various parameters of the method, and show that a fine parameter tuning can lead to significant improvements and stateoftheart performance for GRN inference. Finally, we show that TIGRESS performs well when TFs are not known in advance, i.e. it can predict edge directionality. Overall this study confirms the potential of stateoftheart feature selection techniques for GRN inference.
Methods
Problem formulation
We consider a set of p genes $\mathcal{G}=[1,p]$, including a subset $\mathcal{T}\subset [1,p]$ of transcription factors, among which we wish to discover direct interactions of the form (t,g) for $t\in \mathcal{T}$ and $g\in \mathcal{G}$. We do not try to infer selfregulation, meaning that for each target gene $g\in \mathcal{G}$ we define the set of possible regulators as ${\mathcal{T}}_{g}=\mathcal{T}\setminus \left\{g\right\}$ if $g\in \mathcal{T}$ is itself a transcription factor, and ${\mathcal{T}}_{g}=\mathcal{T}$ otherwise. The set of all candidate regulations is therefore $\mathcal{E}=\left\{(t,g),g\in \mathcal{G},t\in {\mathcal{T}}_{g}\right\}$, and the GRN inference problem is to identify a subset of true regulations among $\mathcal{E}$.
For that purpose, we assume we have gene expression measurements for all genes $\mathcal{G}$ in n experimental conditions. Although the nature of the experiments may vary and typically include knockdown or knockout experiments and even replicates, for simplicity we do not exploit this information and only consider the n×p data matrix of expression levels X as input for the GRN inference problem. Each row of X corresponds to an experiment and each column to a gene. We assume that the expression data have been preprocessed for quality control and missing values imputation.
In order to infer the regulatory network from the expression data X, we compute a score $s:\mathcal{E}\to \mathbb{R}$ to assess the evidence that each candidate regulation is true, and then predict as true regulation the pairs $(t,g)\in \mathcal{E}$ for which the evidence s(t,g) is larger than a threshold δ. We let δ as a usercontrolled parameter, where larger δvalues correspond to less predicted regulations, and only focus on designing a significance score s(t,g) that leads to “good” prediction for some values of δ. In other words, we only focus on finding a good ranking of the candidate regulations $\mathcal{E}$, by decreasing score, such that true regulations tend to be at the top of the list; we let the user control the level of false positive and false negative predictions he can accept. Note that such a ranking is the standard prediction format of the DREAM challenge.
GRN inference with feature selection methods
Many popular methods for GRN inference are based on such a score. For example, the correlation or mutual information between the expression levels of t and g along the different experiments is a popular way to score candidate regulations [12–14]. A drawback of such direct approaches is that it is then difficult to separate direct from indirect regulations. For example, if t_{1}regulates t_{2}which itself regulates g, then the correlation or mutual information between t_{1} and g is likely to be large, although (t_{1}g) is not a direct regulation. Similarly, if t_{1}regulates both t_{2}and g, then t_{2}and g will probably be very correlated, even if there is no direct regulation between them. In order to overcome this problem, a possible strategy is to postprocess the predicted regulations and try to remove regulations likely to be indirect because they are already explained by other regulations [13]. Another strategy is, given a target gene $g\in \mathcal{G}$, to jointly estimate the scores s(t g) for all candidate regulators $t\in {\mathcal{T}}_{g}$ simultaneously, with a method able to capture the fact that a large score for a candidate regulation (t g) is not needed if the apparent correlation between t and g is already explained by other, more likely regulations.
Mathematically, the latter strategy is closely related to the problem of feature selection in statistics, as already observed and exploited by several authors [20, 23]. More specifically, for each target gene $g\in \mathcal{G}$, we consider the regression problem where we wish to predict the expression level of g from the expression level of its candidate regulators $t\in {\mathcal{T}}_{g}$:
where X_{ i }represents the expression level of the ith gene across different experiments (modeled as a random variable), ${X}_{{\mathcal{T}}_{g}}=\left\{{X}_{t}\phantom{\rule{0.3em}{0ex}},\phantom{\rule{0.3em}{0ex}}t\in {\mathcal{T}}_{g}\right\}$ is the set of expression levels of the candidate transcription factors for gene g, and ε is some noise. Any linear or nonlinear statistical method for regression can potentially be used to infer f_{ g }from the observed expression data. However, we are not directly interested in the regression function f_{ g }, but instead in the identification of a small set of transcription factors which are sufficient to provide a good model for X_{ g }. We therefore need a score s_{ g }(t) for each candidate transcription factor $t\in {\mathcal{T}}_{g}$ to assess how likely it is to be involved in the regression model f_{ g }. For example, if we model f_{ g }as a linear function
then the score s_{ g }(t) should typically assess the probability that β_{t,g}is nonzero [23]. More general models are possible, for example [20] model f_{ g } with a random forest [29] and score a predictor s_{ g }(t)with a variable importance measure specific to this model. Once a score s_{ g }(t) is chosen to assess the significance of each transcription factor in the targetgenespecific regression model (1), we can combine them across all target genes by defining the score of a candidate regulation $(t,g)\in \mathcal{E}$ as s(t g)=s_{ g }(t), and rank all candidate regulations by decreasing score for GRN inference.
Feature selection with LARS and stability selection
We now propose a new scoring function s_{ g }(t) to assess the significance of a transcription factor $t\in {\mathcal{T}}_{g}$ in the regression model (1). Our starting point to define the scoring function is the LARS method for feature selection in regression [24]. LARS models the regression function (1) linearly, i.e. it models the expression of a target gene as a linear combination of the expression of its transcription factors, as in (2). Starting from a constant model where no TF is used, it iteratively adds TF in the model to refine the prediction of X_{ g }. Contrary to classical forward stepwise feature selection [30, 31], LARS does not fully reoptimize the fitted model when a new TF is added to the model, but only refines it partially. This results in a statistically sound procedure for feature selection, akin to forward stagewise linear regression and the Lasso [27, 31], and a very efficient computational procedure. In practice, after L steps of the LARS iteration, we obtain a ranked list of L TF selected for their ability to predict the expression of the target gene of interest. Efficient implementations of LARS exist in various programming languages including R (lars package, [24]) and MATLAB (SPAMS toolbox, [32]). Since the selection of TF is iterative, LARS has the potential to disregard indirect regulations.
The direct use of LARS to score candidate regulations has, however, two shortcomings. First, LARS can be very sensitive and unstable in terms of selected features when there exist high correlations between different explanatory variables. Second, it only provides a ranking of the TF, for each TG of interest, but does not provide a score s_{ g }(t) to quantify the evidence that a TF t regulates a target gene g. Since we want to aggregate the predicted regulations across all target genes to obtain a global ranking of all candidate regulations, we need such a score.
To overcome both issues, we do not directly score candidate regulations with the LARS, but instead perform a procedure known as stability selection[26] on top of LARS. The general idea of stability selection is to run a feature selection method many times on randomly perturbed data, and score each feature by the number of times it was selected. It was shown that stability selection can reduce the sensitivity of LARS and Lasso to correlated features, and improve their ability to select correct features [25, 26]. In addition, it provides a score for each feature, which can then be aggregated over different regression problems, i.e. different target genes in our case. More precisely, to score the candidate target genes $t\in {\mathcal{T}}_{g}$ of a given target gene g using LARS with stability selection, we fix a (large) number of iterations R, and repeat R/2times the following iterations: we randomly split the experiments into two halves of equal or approximately equal size, we multiply the expression levels of the candidate transcription factors in ${\mathcal{T}}_{g}$ on each microarray by a random number uniformly sampled on the interval α,1] for some 0≤α≤1, and we run the LARS method for L>0steps on the two resulting reduced and reweighed expression matrices. We therefore perform a total of R LARS runs on randomly modified expression matrices. For each run, the result of LARS after L steps is a ranked list of L TF. After the R runs, we record for each $g\in \mathcal{G}$, $t\in {\mathcal{T}}_{g}$ and l∈[1,Lthe frequency F(g t l) with which the TF t was selected by the LARS in the top l features to predict the expression of gene g. We thus obtain a final score between 0 and 1, 1 meaning that t is always selected by LARS in the top l features to predict the expression level of g, and 0 that is is never selected. Figure 1 displays graphically these frequencies, for a given gene g fixed, all candidate TF in ${\mathcal{T}}_{g}$, and l=1,…,15. When l increases, the frequency F(g t l) for fixed g and t is nondecreasing because the LARS method selects increasing sets of TF at each step. In addition, since the total number of TF selected after l LARS steps is always equal to l, taking the average over the R LARS runs leads to the equality $\sum _{t\in {\mathcal{T}}_{g}}F(g,t,l)=l$, for any gene g and LARS step l.
Once the frequency table F(g t l) is computed for l=1,…,L, we need to convert it into a unique score s(t g) for each candidate pair (t g). The original stability selection score [25, 26] is simply defined as the frequency of selection in the top L variables, i.e.,
As suggested by Figure 1, this score may be very sensitive to the choice of L. In particular, if L is too small, many TF may have zero score (because there are never selected in the top L TFs), but when L is too large, several TF may have the same score 1 because they are always selected in the top L TFs. To alleviate this possible difficulty, we propose as an alternative score to measure the area under each curve up to L steps, i.e. to consider the following area score:
It is worth noting that for a given target gene g, the sum of the scores over the potential transcription factors does not depend on g. Indeed, for any fixed g, there are exactly L TF selected in the first L LARS steps on any randomly modified expression matrix, which implies that the frequencies of selection also sum to L:
Moreover, the area score is also normalized as follows:
This shows that the scores output by TIGRESS are naturally normalized per target gene, and we therefore do not consider further normalization before aggregating all scores together across target genes.
The difference between s_{ original }(t,g)and s_{ area }(t,g)becomes clear if we consider the rank of t in the list produced by LARS in one run as a random variable H_{ t } (with H_{ t }=1meaning that t is ranked first by LARS). F(g,t,l) is then, by definition, the empirical probability P(H_{ t }≤l) that H_{ t }is not larger than l. The original score has therefore an obvious interpretation as P(H_{ t }≤L), which we can rewrite as:
Interestingly a small computation shows that the area score has a similar probabilistic interpretation:
with
In other words, both the original and the area scores can be expressed as E[ϕ(H_{ t })], although with a different function ϕ. While the original score only assesses how often a feature ranks in the top L, the area score additionally takes into account the value of the rank, with features more rewarded if they are not only in the top L but also frequently with a small rank among the top L. Since s_{ area }integrates the frequency information over the full LARS path up to L steps, it should be less sensitive than s_{ original }to the precise choice of L, and should allow to investigate larger values of L without saturation effects when several curves hit the maximal frequency of 1. We note that other scores of the form E[ϕ(H_{ t })] for nonincreasing function ϕcould be investigated as well.
Parameters of TIGRESS
In summary, the full procedure for scoring all candidate edges in $\mathcal{E}$, which we call TIGRESS, splits the GRN inference problem into p independent regression problems taking each target gene $g\in \mathcal{G}$ in turn, and scores each candidate regulation (t g) for a candidate TF $t\in {\mathcal{T}}_{g}$ with the original (3) or area (4) stability score applied to LARS feature selection. In addition to the choice of the scoring method (original or area), the parameters of TIGRESS are (i) the number of runs R performed in stability selection to compute the scores, (ii) the number of LARS steps L, and (iii) the parameter α∈[0,1]which controls the random reweighting of each expression array in each stability selection run. Apart from R that should be taken as large as possible to ensure that frequencies are correctly estimated, and is only limited by the computational time we can afford to run TIGRESS, the influence of αand L on the final performance of the method are not obvious. Taking α=1means that no weight randomization is performed on the different expression arrays, while α=0 leads to maximal randomization. [26] advocate that a value between 0.2and 0.8 is often a good choice. Regarding the choice of L, [26] mentions that it has usually little influence on the result, but as discussed above, the choice of a good range of values may not be trivial in particular for the original score. We investigate below in detail how the performance of TIGRESS depends on the scoring method and on these parameters R, α and L.
Performance evaluation
We experimentally compare TIGRESS to several other GRN inference methods. We use the MATLAB implementations of CLR [14] and GENIE3 [20]. We run ARACNE [13] using the R package minet. We keep default parameter values for each of these methods. Results borrowed from the DREAM5 challenge [28] were directly obtained by each participating team.
Given a gene expression data matrix, each GRN inference method outputs a ranked list of putative regulatory interactions. Taking only the top K predictions in this list, we can compare them to known regulations to assess the number of true positives (TP, the number of known regulations in the top K predictions), false positives (FP, the number of predicted regulations in the top K which are not known regulations), false negatives (FN, the number of known interactions which are not in the top K predictions) and true negatives (TN, the number of pairs not in the top K predictions which are not known regulations). We then compute classical statistics to summarize these numbers for a given K, including precision (TP/(TP + FP)), recall (TP/(TP + FN)), and fallout (FP/(FP + TN)). We assess globally how these statistics vary with K by plotting the receiver operating characteristic (ROC) curve (recall as a function of fallout) and the precisionrecall curve (precision as a function of recall), and computing the area under these curves (respectively AUROC and AUPR) normalized between 0 and 1.
For the datasets of DREAM5, we further compute a Pvalue for the AUROC and AUPR scores, based on all DREAM5 participants’ predictions. This method, which was used by the DREAM5 organizers to rank the teams, involves randomly drawing edges from the teams’ prediction lists and computing the probabilities of obtaining an equal or larger AUPR (resp. AUROC) by chance. More precisely, random lists are constructed as follows: for each row of the predicted list, an edge at the same position is drawn at random from all predictions. For an ensemble of such random lists, the areas under the curves are computed, allowing to estimate a random distribution. Pvalues were obtained by extrapolating the resulting histogram. We refer to [28] for more details on this scoring scheme. Finally, we compute a score for a GRN inference method by integrating the AUROC and AUPR Pvalues as follows:
Data
We evaluate the performance of TIGRESS and other GRN inference methods on nine benchmark datasets, each consisting of a compendium of gene expression data, optionally a list of known TF, and a gold standard set of verified interactions which we ideally would like to recover from the expression data only. Expression data are either simulated or experimentally measured under a wide range of genetic, drug and environmental perturbations. Table 1 summarizes the statistics of these nine networks.
The first three benchmarks are taken from the DREAM5 challenge [28]. Network 1 is a simulated dataset. Its topology and dynamics were modeled according to known GRN, and the expression data were simulated using the GeneNetWeaver software [33]. We refer the interested reader to [22, 34] for more information on this network. The second and third benchmarks are Network 3 and Network 4 of the DREAM5 competition, corresponding respectively to real expression data for E. coli and S. cerevisiae. Note that we do not use in our experiments Network 2 of DREAM5, because no verified TFTG interaction is provided for this dataset consisting in expression data for S. aureus.
Additionally, we run experiments on the E. coli dataset from [14], which has been widely used as a benchmark in GRN inference literature. The expression data was downloaded from the Many Microbe Microarrays (M^{3D}) database [35] (version 4 build 6). It consists in 907 experiments and 4297 genes. We obtained the gold standard data from RegulonDB [36] (version 7.2, May 6th, 2011) that contains 3812 verified interactions among 1525 of the genes present in the microarrays experiments.
Finally, we borrowed the five DREAM4 [22] size 100 multifactorial networks [34] for which the TFs are not known in advance in order to assess TIGRESS’ ability to predict directionality.
As a preprocessing step, we simply meancenter and scale to unit variance the expression levels of each gene within each compendium.
Results
DREAM5 challenge results
In 2010 we participated to the DREAM5 Network Inference Challenge, an open competition to assess the performance of GRN methods [28]. Participants were asked to submit a ranked list of predicted interactions from four matrices of gene expression data. At the time of submission, no further information was available to participants (besides the list of TF), in particular the “true” network of verified interactions for each dataset was not given. After submissions were closed, the organizers of the challenge announced that one network (Network 1) was a simulated network with simulated expression data, while the other expression datasets were real expression data collected for E. coli (Network 3) and S. cerevisiae (Network 4), respectively. Teams were ranked for each network by decreasing score (5), and an overall score was computed summarizing the networkspecific pvalues [28].
We submitted predictions for all networks with a version of TIGRESS that we could not optimize since the benchmarks were blinded at the time of the challenge. We refer to it as Naive TIGRESS below. Naive TIGRESS is the variant of TIGRESS which scores candidate interactions with the original score (3) and uses the arbitrarily fixed parameters α=0.2, L=5, R_{1}=4,000, R_{3}=R_{4}=1,000, where R_{ i } refers to the number of runs for network i. The number of runs were simply set to ensure that TIGRESS would finish within 2 days on a singlecore laptop computer. R_{1}is larger than R_{3}and R_{4}because the size of network 1 is smaller than that of networks 3 and 4, implying that each TIGRESS run is faster. The choice α=0.2followed previous suggestions for the use of stability selection [26], while the choice L=5 roughly corresponded to the largest value for which no TFTG pair had a score of 1.
Naive TIGRESS was among the top GRN prediction methods at DREAM5, ranking second among 29 participating teams in the in silico network challenge, and third overall. Table 2 summarizes the results of the first three teams as well as two stateoftheart methods in average overall score.
The winning method, both in silico and overall, was the GENIE3 method of [20]. GENIE3 already won the DREAM4 challenge, confirming its overall stateoftheart performance. It had particularly strong performance on the in silico network, and more modest performance on both in vivo networks. The ANOVAbased method of [37] ranked second overall, with particularly strong performance on the E. coli network. Naive TIGRESS ranked third overall, with particularly strong performance on the in silico network, improving over GENIE3 in terms of AUPR.
Interestingly, GENIE3 and TIGRESS follow a similar formulation of GRN inference as a collection of feature selection problems for each target gene, and use similar randomizationbased techniques to score the evidence of a candidate TFTG interaction. The main difference between the two methods is that GENIE3 aggregates the features selected by decision trees, while TIGRESS aggregates the features selected by LARS. The overall good results obtained by both methods suggest that this formalism is particularly relevant for GRN inference.
Influence of TIGRESS parameters
In this section, we provide more details about the influence of the various parameters of TIGRESS on its performance, taking DREAM5 in silico network as benchmark dataset. Obviously the improvements we report below would require confirmation on new datasets not used to optimize the parameters, but they shed light on the further potential of TIGRESS and similar regressionbased method when parameters are precisely tuned.
Starting from the parameters used in Naive TIGRESS (R=4,000, α=0.2and L=5, original score), we assess the influence of the different parameters by systematically testing the following combinations:

original (3) or area (4) scoring method;

randomization parameter α∈{0,0.1…,1};

length of the LARS path L∈{1,2…20};

number of randomization runs R∈{1,000;4,000;10,000}.
Figure 2 summarizes the score (5) obtained by each combination of parameters on Network 1.
A first observation is that the area scoring method consistently outperforms the original scoring method, for any choice of αand L. This suggests that, by default, the newly proposed area score should be preferred to the classical original score. We also note that the performance of the area score is less sensitive to the value of α or L than that of the original score. For example, any value of αbetween 0.2and 0.8, and any L less than 10 leads to a score of at least 90 for the area score, but it can go down to 60 for the original score. This is a second argument in favor of the area scoring setting: as it is not very sensitive to the choice of the parameters, one may practically more easily tune it for optimal performance. On the contrary, the window of (α L) values leading to the best performance is more narrow with the original scoring method, and therefore more difficult to find a priori. The recommendation of [26] to choose αin the range [0.2,0.8]is clearly not precise enough for GRN inference. The best overall performance is obtained with (α=0.4,L=2) in both scoring settings.
Regarding the relationship between αand L, we observe in Figure 2 a slight positive correlation for the optimal L as a function of α, particularly for the area score. For example, for R=10^{4}, L=2 is optimal for α≤0.4, but L≥4 is optimal for α≥0.8. The effect is even more pronounced for R=4,000. This can be explained by the fact that when αincreases, we decrease the variations in the the different runs of LARS and therefore reduce the diversity of features selected; increasing the number of LARS is a way to compensate this effect by increasing the number of features selected at each run. Another way to observe the need to ensure a sufficient diversity is to observe how the best parameters L and αvary as a function of R (Figure 3). It appears clearly that the optimal number of steps L^{∗} decreases when the number of resampling runs increases and stabilizes at 2. This is not a surprising result. Indeed, when more resampling is performed, the chance of selecting a given feature increases. The number N of non zero scores subsequently increases and it thus becomes unnecessary to look further in the regularization path. On the other hand, the value of α^{∗} lies steadily between 0.3and 0.5, suggesting that the adjustment to the number of bootstraps can mostly be made through the choice of L.
Furthermore, we unsurprisingly observe that increasing the number R of resampling runs leads to better performances. On Figure 4, we show the score as a function of R with L=2 and α=0.4. We clearly see that, for both scoring methods, increasing the number of runs is beneficial. The performance seems to reach an asymptote only when R becomes larger than 5,000.
Finally, we were interested in the number of TFs selected per gene. Figures 5 and 6 show how the distribution of this number changes with respect to the total number of predictions for L=2 and L=20respectively. We observe a lower variance and a larger median when L is larger, which suggests that choosing a small value for L leads to predicting more variable numbers of interactions per TG whereas a large value will force all TGs to be linked to a similar and higher number of TFs. This observation sheds some light on the choice of L in general, when assumptions can be made on the topology of the network to predict.
Comparison with other methods
Figure 7 depicts both the ROC and the Precision/Recall curves for several methods on Network 1. Table 2 summarizes these performances in terms of AUPR, AUROC and related pvalues as well as the score (5). Here, TIGRESS was run with α=0.4, L=2 and R=8,000 which corresponds to the best performance of the algorithm, as investigated in the previous section.
TIGRESS, as tuned optimally on this network, outperforms all methods in terms of AUPR and all methods but GENIE3 in terms of AUROC. Moreover, the shape of the Precision/Recall curve suggests that the top of the prediction list provided by TIGRESS contains more true edges than other methods. The ROC curve, on the other hand, focuses on the entire list of results. Therefore, we would argue that TIGRESS can be more reliable than GENIE in its first predictions but contains overall more errors when we go further in the list.
These results suggest that TIGRESS has the potential to compare with stateoftheart methods and confirm the importance of correct parameter tuning.
In vivo networks results
Since Naive TIGRESS did not perform very well on the in vivo networks at the DREAM5 competition (Table 2), we now test on these networks TIGRESS with the best parameters selected on the in silico (area score, α=0.4, L=2 and R=10,000). Table 2 also shows the values of AUPR, AUROC, related pvalues and score for DREAM5 networks 3 and 4 reached by TIGRESS and ROC and P/R curves were drawn on Figures 8 and 9.
The results on these two networks are overall disappointing: TIGRESS does not do better than Naive TIGRESS. In fact, both sets of results are very weak. Without attempting to reoptimize all parameters for each network, one may wonder whether the parameters chosen using the in silico network are optimal for the in vivo networks. As a partial answer, Figure 10 shows the behavior of the score with respect to L for Networks 3 and 4. Interestingly, it seems that a much larger L is preferable in this case, suggesting that one may have to adapt the parameters to the size of the network in terms of number of transcription factors. Indeed, networks 3 and 4 contain respectively 334 and 333 transcription factors, making them much larger than the in silico and the E. coli networks (195 and 180 TFs respectively), for which a small L leads to a better performance. Choosing L=100 for DREAM5 in vivo networks yields much better results. As a matter of fact, TIGRESS obtains the best results on Network 4 with this value of L and doubles its performance on Network 3.
On Figure 11 we compare Precision/Recall and ROC curves obtained with TIGRESS with several other algorithms on the E. coli network from^{1}[14]. Table 3 compares the areas under the curves. TIGRESS is comparable to CLR, while GENIE3 outperforms other methods. However the overall performance of all methods remains disappointing.
Analysis of errors on E. coli
To understand further the advantages and limitations of TIGRESS, we analyze the type of errors it typically makes taking the E. coli dataset as example. We analyze FP, i.e. cases where TIGRESS predicts an interaction that does not appear in the gold standard GRN.
We focus in particular on quantifying how far a wrongly predicted interaction is from a true one, and introduce for that purpose the notion of distance between two genes as the shortest path distance between them on the gold standard GRN, forgetting about the direction of edges. For two genes G1 and G2, we call G1G2 a distancex link if the shortest path between G1 and G2 on the true network has length x. Figure 12 shows the distribution of these distances for spuriously discovered edges over the gold standard network, i.e. the actual proportion of distancex links, with x∈{1,2,3,4,>4}. We write ${\widehat{p}}_{x}$ the proportion of spurious TFTG interactions with distance x.
Figure 13 depicts the distribution of distancex proportions among the spuriously detected edges, as a function of the number of predicted edges. Dotted lines represent the 95% confidence interval around the exact distribution${\left({\widehat{p}}_{x}\right)}_{x}$. For a given number r of spuriously predicted edges, this interval is computed as
where${q}_{a}\left({\widehat{p}}_{x}\right)$ represents the quantile of order a of a hypergeometric distribution$\mathcal{H}({N}_{S},{\widehat{p}}_{x}{N}_{S},r)$ and N_{ S } is the total number of spuriously predicted edges.
We observe that most of the recovered false positives appear as distance2 edges in a significantly higher proportion than${\widehat{p}}_{2}$ whereas significantly less distance>4edges are discovered. These results strongly suggest that most of TIGRESS errors  especially at the top of the list  are indeed sensible guesses, where the two nodes, spuriously discovered with a parent/child relationship are actually separated by only one other node. In Figure 14, we detail the three possible patterns observable in this situation.
Figure 15 focuses on distance2 errors. Note that some edges show more than one pattern, e.g. the first spurious edges are both siblings and couples. It appears that most of them are siblings and can thus be interpreted as spurious feedforward loops. We believe that this can be explained by three main reasons: i) the discovered edges actually exist but have not been experimentally validated yet; ii) there is more of a linear relationship between siblings than between parent and child; iii) some nodes have very correlated expression levels, making it difficult for TIGRESS to tell between the parent and the child.
Directionality prediction : case study on DREAM4 networks
In order to check whether TIGRESS can predict edge directions, we additionally ran it on the five size 100 multifactorial DREAM4 networks, for which the TFs are not known. The five datasets contain 100 samples and 100 genes. We observe that TIGRESS can indeed perform well in this setting. Table 4 shows the results with the default parameter setting (L=2, α=0.4, R=10,000) compared with those of GENIE3, that won the DREAM4 network inference challenge. Without further optimization of the parameters on these networks, TIGRESS achieves a better overall performance than GENIE3.
Furthermore, we ran the complete analysis of the parameters on these networks, to check whether the optimal parameters change when TFs are unknown. Figure 16 shows the overall performance, that is the average performance on all five networks, as a function of α and L. Given the small size of the networks, it is not surprising that the optimalL is equal to 1. It also seems that the optimal value forα lies in between 0 and 0.1, corresponding to a strong randomization.
Computational Complexity
The complexity of runningL LARS steps on a regression problem withq covariates and n samples isO(nqL + L^{3})[24]. In our case,q is the number of TF andn is the number of expression arrays, which we divide by two during the resampling step, and we pay this complexity for each TG and each resampling. Multiplying byp TG andR resampling runs, we therefore get a total complexity of orderO(pR(n /2 qL +L^{3})), which boils down toO(pRnqL/2) in the situation whereL is smaller thann/2andq.
Table 5 compares the running time of TIGRESS, GENIE3, ARACNE, ANOVerence and CLR on Network 1. We define runningunits units for GENIE3 and TIGRESS that correspond to one element of computation, i.e. the equivalent ofK=T=p=q=n=1 for GENIE3 andL=R=p=q=n=1 for TIGRESS. As stated in [20], the complexity of GENIE3 is of the order ofO(pTKn log(n)). It is then a matter of multiplication to get the approximate running time for a given dataset and a given set of parameters. The total running time for these two methods is computed using default parameters, that is$\mathit{\text{GENIE}}3(T=1,000;K=\sqrt{q})$ and TIGRESS(L=2;R=10,000). All algorithms were run on a 12GB RAM Intel X5472 3.00GHz computer.
Discussion and conclusions
In this paper, we presented TIGRESS, a new method for GRN inference. TIGRESS expresses the GRN inference problem as a feature selection problem, and solves it with the popular LARS feature selection method combined with stability selection. It ranked in the top 3 GRN inference methods at the 2010 DREAM5 challenge, without any parameter tuning. We clarified in this paper the influence of each parameter, and showed that further improvement may result from finer parameter tuning.
We proposed in particular a new scoring method for stability selection, based on the area under the stability curve. It differs from the original formulation of [26] which does not take into account the full distribution of ranks of a feature in the randomized feature selection procedure. Comparing the two, we observed that the new area scoring technique yields better results and is less sensitive to the values of the parameters: practically all values of,e.g., the randomization parameter α yield the same performance. Similarly, the choice of the number L of LARS steps to run seems to have much less impact on the performance in this new setting. As we showed, the original and area scores for a featuret can be both expressed in a common formalism asE ϕ(H)] for different functionsϕ, where H_{ t }is the rank of feature t as selected by the randomized LARS. It could be interesting to systematically investigate variants of these scores with more general nonincreasing functionsϕ, not only for GRN inference but also more generally as a generic feature selection procedure.
Comparing TIGRESS  as tuned optimally  to stateoftheart algorithms on thein silico network, we observed that it achieves a similar performance to that of GENIE3 [20], the best performer at the DREAM5 challenge. However, TIGRESS does not do as good as this algorithm onin vivo networks. GENIE3 is also an ensemble algorithm but differs from TIGRESS in that it uses a nonlinear treebased method for feature selection, while TIGRESS uses LARS. The difference in performance could be explained by the fact that the linear relationship between TGs and TFs assumed by TIGRESS is farfetched given the obvious complexity of the problem.
A further analysis of our results on the E. coli network from [14] showed that many spuriously detected edges follow the same pattern: TIGRESS discovers edges when in reality the two nodes aresiblings, and thus tends to wrongly predict feedforward loops. This result suggests many directions for future work. Among them, we believe that operons,i.e. groups of TGs regulated together could be part of the problem. Moreover, it could be that there is more of a linear relationship between siblings than between parent and child, as TFs are known to be operating asswitches, i.e. it is only after a certain amount change in expression of the TF that related TGs are affected. However, it is worth noting thatin vivo networks gold standards may not be complete. Therefore, the hypothesis that TIGRESS is actually correct when predicting these loops cannot be discarded.
TIGRESS depends on four parameters: the scoring method, the numberR of resampling runs, the randomization factorαand the number of LARS stepsL. We showed in this paper that changing the value of these parameters can greatly affect the performance and provided guidelines to choose them. It is worth noting, though, that other modifications can be imagined. In particular, one may wonder about the influence of the resampling parameters (with or without replacement, proportion of samples to be resampled). These questions will be tackled in future work.
While it seems indeed more realistic not to restrict underlying models to linear ones, it is fair to say that no method performs very well in absolute values onin vivo networks. For example, performances on theE. coli network seem to level out at some 64% AUROC and 8% AUPR which cannot be considered satisfying. This suggests that while regressionbased procedures such as TIGRESS or GENIE3 are stateoftheart for GRN inference, their performances seem to hit a limit which probably cannot be outdistanced without some changes. It is worth noting that, as argued in [28], combining these methods together leads to improvement, as different sets of interactions are discovered by each method. Another way to overcome these limits may be a change in the global approach such as adding some supervision in the learning process as,e.g., investigated in [38].
Endnotes
^{1} Unfortunately, we were not able to run ANOVerence on this particular dataset due to heavy formatting requirements of the data by the algorithm that we did not have the necessary information to perform.
References
 1.
Arkin A, Shen P, Ross J: A test case of correlation metric construction of a reaction pathway from measurements. Science. 1997, 277 (5330): 12751279. 10.1126/science.277.5330.1275. [http://www.sciencemag.org/cgi/reprint/277/5330/1275.pdf] 10.1126/science.277.5330.1275
 2.
Liang S, Fuhrman S, Somogyi R: REVEAL, a general reverse engineering algorithm for inference of genetic network architectures. Pac Symp Biocomput. 1998, 3: 1829.
 3.
Chen T, He HL, Church GM: Modeling gene expression with differential equations. Pac Symp Biocomput. 1999, 4: 2940.
 4.
Akutsu T, Miyano S, Kuhara S: Algorithms for identifying Boolean networks and related biological networks based on matrix multiplication and fingerprint function. J Comput Biol. 2000, 7 (34): 331343. 10.1089/106652700750050817.
 5.
Yeung MKS, Tegnér J, Collins JJ: Reverse engineering gene networks using singular value decomposition and robust regression. Proc Natl Acad Sci USA. 2002, 99 (9): 61636168. 10.1073/pnas.092576199. [http://www.pnas.org/content/99/9/6163.abstract] 10.1073/pnas.092576199
 6.
Tegner J, Yeung MKS, Hasty J, Collins JJ: Reverse engineering gene networks: integrating genetic perturbations with dynamical modeling. Proc Natl Acad Sci USA. 2003, 100 (10): 59445949. 10.1073/pnas.0933416100.
 7.
Gardner TS, Bernardo D, Lorenz D, Collins JJ: Inferring genetic networks and identifying compound mode of action via expression profiling. Science. 2003, 301 (5629): 102105. 10.1126/science.1081900.
 8.
Chen KC, Wang TY, Tseng HH, Huang CYF, Kao CY: A stochastic differential equation model for quantifying transcriptional regulatory network in Saccharomyces cerevisiae. Bioinformatics. 2005, 21 (12): 28832890. 10.1093/bioinformatics/bti415.
 9.
Bernardo D, Thompson MJ, Gardner TS, Chobot SE, Eastwood EL, Wojtovich AP, Elliott SJ, Schaus SE, Collins JJ: Chemogenomic profiling on a genomewide scale using reverseengineered gene networks. Nat Biotechnol. 2005, 23 (3): 377383. 10.1038/nbt1075.
 10.
Bansal M, Della Gatta, Bernardo D: Inference of gene regulatory networks and compound mode of action from time course gene expression profiles. Bioinformatics. 2006, 22 (7): 815822. 10.1093/bioinformatics/btl003.
 11.
Zoppoli P, Morganella S, Ceccarelli M: TimeDelayARACNE: Reverse engineering of gene networks from timecourse data by an information theoretic approach. BMC Bioinformatics. 2010, 11: 15410.1186/1471210511154.
 12.
Butte AJ, Tamayo P, Slonim D, Golub TR, Kohane IS: Discovering functional relationships between RNA expression and chemotherapeutic susceptibility using relevance networks. Proc Natl Acad Sci USA. 2000, 97 (22): 1218212186. 10.1073/pnas.220392197.
 13.
Margolin AA, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Dalla Favera R, Califano A: ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular contexts. BMC Bioinformatics. 2006, 7 Suppl 1: S710.1186/147121057S1S7.
 14.
Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins JJ, Gardner TS: Largescale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007, 5: e810.1371/journal.pbio.0050008.
 15.
Rice J, Tu Y, Stolovitzky G: Reconstructing biological networks using conditional correlation analysis. Bioinformatics. 2005, 21 (6): 765773. 10.1093/bioinformatics/bti064.
 16.
Friedman N, Linial M, Nachman I, Pe’er D: Using Bayesian networks to analyze expression data. J Comput Biol. 2000, 7 (34): 601620. 10.1089/106652700750050961.
 17.
Hartemink A, Gifford D, Jaakkola T, Young R: Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks. Proceedings of the Pacific Symposium on Biocomputing 2002. Edited by: Altman RB, Dunker AK, Hunter L, Lauerdale K, Klein TE. 2002, World Scientific, 422433. [http://helixweb.stanford.edu/psb01/abstracts/p422.html]
 18.
Perrin B, 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. 10.1093/bioinformatics/btg1071.
 19.
Friedman N: Inferring cellular networks using probabilistic graphical models. Science. 2004, 303 (5659): 79910.1126/science.1094068.
 20.
HuynhThu VA, Irrthum A, Wehenkel L, Geurts P: Inferring regulatory networks from expression data using treebased methods. PLoS One. 2010, 5 (9): e1277610.1371/journal.pone.0012776.
 21.
Markowetz F, Spang R: Inferring cellular networks  a review. BMC Bioinformatics. 2007, 8 (Suppl 6): S510.1186/147121058S6S5. [http://www.biomedcentral.com/14712105/8/S6/S5] 10.1186/147121058S6S5
 22.
Marbach D, Prill RJ, Schaffter T, Mattiussi C, Floreano D, Stolovitzky G: Revealing strengths and weaknesses of methods for gene network inference. Proc Natl Acad Sci USA. 2010, 107 (14): 62866291. 10.1073/pnas.0913357107.http://www.pnas.org/content/107/14/6286.abstract, 10.1073/pnas.0913357107
 23.
Meinshausen N, Bühlmann P: High dimensional graphs and variable selection with the Lasso. Ann Stat. 2006, 34: 14361462. 10.1214/009053606000000281.
 24.
Efron B, Hastie T, Johnstone I, Tibshirani R: Least angle regression. Ann. Stat. 2004, 32 (2): 407499. 10.1214/009053604000000067.
 25.
Bach FR: Bolasso: model consistent Lasso estimation through the bootstrap. Proceedings of theth international conference on Machine learning Volume 308 of ACM International Conference Proceeding Series. Edited by: Cohen WW, McCallum A, Roweis ST. 2008, ACM, New York, NY, USA, 3340.
 26.
Meinshausen N, Bühlmann P: Stability selection. J R Stat Soc Ser B. 2010, 72 (4): 417473. 10.1111/j.14679868.2010.00740.x.
 27.
Tibshirani R: Regression shrinkage and selection via the lasso. J R Stat Soc Ser B. 1996, 58: 267288.
 28.
Marbach D, Costello J, Küffner R, Vega N, Prill R, Camacho D, Allison K, Kellis M, Collins J, Stolovitzky G, the DREAM5 Consortium: Wisdom of crowds for robust gene network inference. Nat Methods. 2012, 9 (8): 796804. 10.1038/nmeth.2016.
 29.
Breiman L: Random forests. Mach Learn. 2001, 45: 532. 10.1023/A:1010933404324.
 30.
Weisberg S: Applied linear regression. 1981, NewYork, Wiley
 31.
Hastie T, Tibshirani R, Friedman J: The elements of statistical learning: data mining, inference, and prediction. 2001
 32.
Mairal J, Bach F, Ponce J, Sapiro G: Online Learning for Matrix Factorization and Sparse Coding. J Mach Learn Res. 2010, 11: 1960. [http://jmlr.csail.mit.edu/papers/v11/mairal10a.html]
 33.
Schaffter T, Marbach D, Floreano D: GeneNetWeaver: in silico benchmark generation and performance profiling of network inference methods. Bioinformatics. 2011, 27 (16): 22632270. 10.1093/bioinformatics/btr373. [http://bioinformatics.oxfordjournals.org/content/27/16/2263.abstract] 10.1093/bioinformatics/btr373
 34.
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. 10.1089/cmb.2008.09TT. [http://online.liebertpub.com/doi/abs/10.1089/cmb.2008.09TT] 10.1089/cmb.2008.09TT
 35.
Faith J, Driscoll M, Fusaro V, Cosgrove E, Hayete B, Juhn F, Schneider S, Gardner T: Many Microbe Microarrays Database: uniformly normalized Affymetrix compendia with structured experimental metadata. Nucleic Acids Res. 2008, 36 (Database issue): D866—D87010.1093/nar/gkm815.
 36.
GamaCastro S, Salgado H, PeraltaGil M, SantosZavaleta A, MuñizRascado L, SolanoLira H, JimenezJacinto V, Weiss V, GarcíaSotelo JS, LópezFuentes A, PorrónSotelo L, AlquiciraHernández S, MedinaRivera A, MartínezFlores I, AlquiciraHernández K, MartínezAdame R, BonavidesMartínez C, MirandaRíos J, Huerta AM, MendozaVargas A, ColladoTorres L, Taboada B, VegaAlvarado L, Olvera M, Olvera L, Grande R, Morett E, ColladoVides J: 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): D98—D105[http://nar.oxfordjournals.org/content/39/suppl_1/D98.abstract]
 37.
Küffner R, Petri T, Tavakkolkhah P, Windhager L, Zimmer R: Inferring gene regulatory networks by ANOVA. Bioinformatics. 2012, 28 (10): 13761382. 10.1093/bioinformatics/bts143.
 38.
Mordelet F, Vert JP: SIRENE: Supervised inference of regulatory networks. Bioinformatics. 2008, 24 (16): i76—i8210.1093/bioinformatics/btn273.
Acknowledgements
JPV was supported by the French National Reseach Agency (ANR09BLAN005104 and ANR11BSV201302) and the European Research Council (SMACERC280032).
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
All authors conceived the method and drafted the manuscript. ACH and JPV implemented the method and analyzed the results. ACH ran the experiments. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Gene Regulatory Network inference
 Feature selection
 Gene expression data
 LARS
 Stability selection