Skip to main content
  • Methodology article
  • Open access
  • Published:

TIGRESS: Trustful Inference of Gene REgulation using Stability Selection

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 regression-based 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 state-of-the-art performance for GRN inference, in both directed and undirected settings.

Conclusions

TIGRESS reaches state-of-the-art 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 (GP-DREAM, 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 TF-TG interactions has many potential far-reaching 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 TF-TG 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 high-throughput methods, in particular DNA microarrays, to monitor gene expression on a genome-wide 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 TF-TG 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 TF-TG 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 [111]. 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 information-theoretic measures such as mutual information [1214] or take explicit advantage of perturbations in the experiments such as gene knock-downs [15]. The difficulty to separate direct from indirect regulations has been addressed with the formalism of Bayesian networks [1619], or by formulating the GRN inference problem as a feature selection problem [20]. Mutual information-based 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 state-of-the-art on several benchmarks [20, 22]. Feature selection with random forests remains however poorly understood theoretically, and one may wonder how other well-established 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 state-of-the-art 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 state-of-the-art feature selection techniques for GRN inference.

Methods

Problem formulation

We consider a set of p genes G=[1,p], including a subset T[1,p] of transcription factors, among which we wish to discover direct interactions of the form (t,g) for tT and gG. We do not try to infer self-regulation, meaning that for each target gene gG we define the set of possible regulators as T g =T{g} if gT is itself a transcription factor, and T g =T otherwise. The set of all candidate regulations is therefore E= ( t , g ) , g G , t T g , and the GRN inference problem is to identify a subset of true regulations among E.

For that purpose, we assume we have gene expression measurements for all genes G in n experimental conditions. Although the nature of the experiments may vary and typically include knock-down or knock-out 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 pre-processed for quality control and missing values imputation.

In order to infer the regulatory network from the expression data X, we compute a score s:ER to assess the evidence that each candidate regulation is true, and then predict as true regulation the pairs (t,g)E for which the evidence s(t,g) is larger than a threshold δ. We let δ as a user-controlled 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 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 [1214]. A drawback of such direct approaches is that it is then difficult to separate direct from indirect regulations. For example, if t1regulates t2which itself regulates g, then the correlation or mutual information between t1 and g is likely to be large, although (t1g) is not a direct regulation. Similarly, if t1regulates both t2and g, then t2and 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 post-process 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 gG, to jointly estimate the scores s(t g) for all candidate regulators t 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 gG, we consider the regression problem where we wish to predict the expression level of g from the expression level of its candidate regulators t T g :

X g = f g ( X T g )+ε,
(1)

where X i represents the expression level of the i-th gene across different experiments (modeled as a random variable), X T g = X t , t T g 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 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

f g ( X T g )= t T g β t , g X t ,
(2)

then the score s g (t) should typically assess the probability that βt,gis non-zero [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 target-gene-specific regression model (1), we can combine them across all target genes by defining the score of a candidate regulation (t,g)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 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 re-optimize 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 stage-wise 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 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 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 gG, t 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 T g , and l=1,…,15. When l increases, the frequency F(g t l) for fixed g and t is non-decreasing 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 t T g F(g,t,l)=l, for any gene g and LARS step l.

Figure 1
figure 1

Stability selection. Illustration of the stability selection frequency F(g,t,L) for a fixed target gene g. Each curve represents a TF t T g , and the horizontal axis represents the number L of LARS steps. F(g,t,L) is the frequency with which t is selected in the first L LARS steps to predict g, when the expression matrix is randomly perturbed by selecting only a limited number of experiments and randomly weighting each expression array. For example, the TF corresponding to the highest curve was selected 57% of the time at the first LARS step, and 81% of the time in the first two LARS steps.

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

s original (t,g)=F(g,t,L).
(3)

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:

s area (t,g)= 1 L l = 1 L F(g,t,l).
(4)

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:

t s original ( t , g ) = t F ( g , t , L ) = L .

Moreover, the area score is also normalized as follows:

t s area ( t , g ) = t 1 L l = 1 L F ( g , t , l ) = 1 L l = 1 L l = L + 1 2 .

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:

s original ( t , g ) = E ϕ original ( H t ) with ϕ original ( h ) = 1 if h L , 0 otherwise.

Interestingly a small computation shows that the area score has a similar probabilistic interpretation:

s area ( t , g ) = l = 1 L F ( g , t , l ) = l = 1 L P ( H t l ) = l = 1 L h = 1 l P ( H t = h ) = h = 1 L ( L + 1 h ) P ( H t = h ) = E ϕ area ( H t ) ,

with

ϕ area ( h ) = L + 1 h if h L , 0 otherwise.

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 non-increasing function ϕcould be investigated as well.

Parameters of TIGRESS

In summary, the full procedure for scoring all candidate edges in E, which we call TIGRESS, splits the GRN inference problem into p independent regression problems taking each target gene gG in turn, and scores each candidate regulation (t g) for a candidate TF t 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 re-weighting 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 fall-out (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 fall-out) and the precision-recall 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 P-value 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. P-values 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 P-values as follows:

overall score= 1 2 log 10 ( P AUPR P AUROC ).
(5)

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.

Table 1 Datasets

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 TF-TG 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 (M3D) 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 pre-processing step, we simply mean-center 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 network-specific p-values [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, R1=4,000, R3=R4=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 single-core laptop computer. R1is larger than R3and R4because 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 TF-TG 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 state-of-the-art methods in average overall score.

Table 2 DREAM5 networks results

The winning method, both in silico and overall, was the GENIE3 method of [20]. GENIE3 already won the DREAM4 challenge, confirming its overall state-of-the-art performance. It had particularly strong performance on the in silico network, and more modest performance on both in vivo networks. The ANOVA-based 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 randomization-based techniques to score the evidence of a candidate TF-TG 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 regression-based 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.

Figure 2
figure 2

Score for network 1. Top plots show the score for R=4,000and bottom plots depict the case R=10,000for both the area (left) and the original (right) scoring settings, as a function of αand L.

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=104, 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.

Figure 3
figure 3

Optimal values of the parameters. Optimal values of parameters L, αand N with respect to the number of resampling runs.

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.

Figure 4
figure 4

Impact of the number of resampling runs. Score as a function of R. In both scoring settings, αand L were set to 0.4and 2, respectively.

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.

Figure 5
figure 5

Distribution of the number of TFs selected per gene for L=2. Histograms of the number of TFs selected per gene with respect to the total number of predictions when L=2.

Figure 6
figure 6

Distribution of the number of TFs selected per gene for L=20. Histograms of the number of TFs selected per gene with respect to the total number of predictions when L=20.

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 p-values 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.

Figure 7
figure 7

Performance on network 1. ROC (left) and Precision/Recall (right) curves for several methods on Network 1.

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 state-of-the-art 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 p-values and score for DREAM5 networks 3 and 4 reached by TIGRESS and ROC and P/R curves were drawn on Figures 8 and 9.

Figure 8
figure 8

Performance on DREAM5 network 3. ROC (Left) and Precision/Recall (Right) curves for several methods on DREAM5 network 3.

Figure 9
figure 9

Performance on DREAM5 network 4. ROC (Left) and Precision/Recall (Right) curves for several methods on DREAM5 network 4.

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 re-optimize 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.

Figure 10
figure 10

In vivo networks results. Score with respect to L for DREAM5 networks 3 and 4 and E. coli network (α=0.4, R=10,000).

On Figure 11 we compare Precision/Recall and ROC curves obtained with TIGRESS with several other algorithms on the E. coli network from1[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.

Figure 11
figure 11

Performance on the E. coli network. ROC (Left) and Precision/Recall (Right) curves for several methods on the E. coli dataset.

Table 3 E. coli network results

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 G1-G2 a distance-x 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 distance-x links, with x{1,2,3,4,>4}. We write p ̂ x the proportion of spurious TF-TG interactions with distance x.

Figure 12
figure 12

Spurious edges shortest path distribution. Exact distribution of the shortest path between spuriously predicted TF-TG couples.

Figure 13 depicts the distribution of distance-x 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 ( p ̂ x ) x . For a given number r of spuriously predicted edges, this interval is computed as

q 0 . 025 ( p ̂ x ) r ; q 0 . 975 ( p ̂ x ) r ,

where q a ( p ̂ x ) represents the quantile of order a of a hypergeometric distributionH( N S , p ̂ x N S ,r) and N S is the total number of spuriously predicted edges.

Figure 13
figure 13

Distribution of the shortest path with respect to the number of predictions. Distribution of the shortest path length between nodes of spuriously detected edges and 95% confidence interval for the null distribution. These edges are ranked by order of discovery.

We observe that most of the recovered false positives appear as distance-2 edges in a significantly higher proportion than 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 14
figure 14

Distance-2 patterns. The three possible distance-2 patterns: siblings, couple and grandparent/grandchild relationships.

Figure 15 focuses on distance-2 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 feed-forward 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.

Figure 15
figure 15

Distribution of distance-2 errors. Distribution of distance 2 errors with respect to the number of predictions. 95% error bars were computed using the quantiles of a hypergeometric distribution.

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.

Table 4 DREAM4 networks results

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.

Figure 16
figure 16

Results on DREAM4 networks. Overall score on the five multifactorial size 100 DREAM4 networks, as a function ofαand L.

Computational Complexity

The complexity of runningL LARS steps on a regression problem withq covariates and n samples isO(nqL + L3)[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 +L3)), 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 isGENIE3(T=1,000;K= q ) and TIGRESS(L=2;R=10,000). All algorithms were run on a 12GB RAM Intel X5472 3.00GHz computer.

Table 5 Runtime

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 non-increasing functionsϕ, not only for GRN inference but also more generally as a generic feature selection procedure.

Comparing TIGRESS - as tuned optimally - to state-of-the-art 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 non-linear tree-based 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 far-fetched 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 feed-forward 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 regression-based procedures such as TIGRESS or GENIE3 are state-of-the-art 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): 1275-1279. 10.1126/science.277.5330.1275. [http://www.sciencemag.org/cgi/reprint/277/5330/1275.pdf] 10.1126/science.277.5330.1275

    Article  CAS  Google Scholar 

  2. Liang S, Fuhrman S, Somogyi R: REVEAL, a general reverse engineering algorithm for inference of genetic network architectures. Pac Symp Biocomput. 1998, 3: 18-29.

    Google Scholar 

  3. Chen T, He HL, Church GM: Modeling gene expression with differential equations. Pac Symp Biocomput. 1999, 4: 29-40.

    Google Scholar 

  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 (3-4): 331-343. 10.1089/106652700750050817.

    Article  CAS  Google Scholar 

  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): 6163-6168. 10.1073/pnas.092576199. [http://www.pnas.org/content/99/9/6163.abstract] 10.1073/pnas.092576199

    Article  CAS  Google Scholar 

  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): 5944-5949. 10.1073/pnas.0933416100.

    Article  CAS  Google Scholar 

  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): 102-105. 10.1126/science.1081900.

    Article  CAS  Google Scholar 

  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): 2883-2890. 10.1093/bioinformatics/bti415.

    Article  CAS  Google Scholar 

  9. Bernardo D, Thompson MJ, Gardner TS, Chobot SE, Eastwood EL, Wojtovich AP, Elliott SJ, Schaus SE, Collins JJ: Chemogenomic profiling on a genome-wide scale using reverse-engineered gene networks. Nat Biotechnol. 2005, 23 (3): 377-383. 10.1038/nbt1075.

    Article  Google Scholar 

  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): 815-822. 10.1093/bioinformatics/btl003.

    Article  CAS  Google Scholar 

  11. Zoppoli P, Morganella S, Ceccarelli M: TimeDelay-ARACNE: Reverse engineering of gene networks from time-course data by an information theoretic approach. BMC Bioinformatics. 2010, 11: 154-10.1186/1471-2105-11-154.

    Article  Google Scholar 

  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): 12182-12186. 10.1073/pnas.220392197.

    Article  CAS  Google Scholar 

  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: S7-10.1186/1471-2105-7-S1-S7.

    Article  Google Scholar 

  14. Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins JJ, Gardner TS: Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007, 5: e8-10.1371/journal.pbio.0050008.

    Article  Google Scholar 

  15. Rice J, Tu Y, Stolovitzky G: Reconstructing biological networks using conditional correlation analysis. Bioinformatics. 2005, 21 (6): 765-773. 10.1093/bioinformatics/bti064.

    Article  CAS  Google Scholar 

  16. Friedman N, Linial M, Nachman I, Pe’er D: Using Bayesian networks to analyze expression data. J Comput Biol. 2000, 7 (3-4): 601-620. 10.1089/106652700750050961.

    Article  CAS  Google Scholar 

  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, 422-433. [http://helix-web.stanford.edu/psb01/abstracts/p422.html]

    Google Scholar 

  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): ii138-ii148. 10.1093/bioinformatics/btg1071.

    Article  Google Scholar 

  19. Friedman N: Inferring cellular networks using probabilistic graphical models. Science. 2004, 303 (5659): 799-10.1126/science.1094068.

    Article  CAS  Google Scholar 

  20. Huynh-Thu VA, Irrthum A, Wehenkel L, Geurts P: Inferring regulatory networks from expression data using tree-based methods. PLoS One. 2010, 5 (9): e12776-10.1371/journal.pone.0012776.

    Article  Google Scholar 

  21. Markowetz F, Spang R: Inferring cellular networks - a review. BMC Bioinformatics. 2007, 8 (Suppl 6): S5-10.1186/1471-2105-8-S6-S5. [http://www.biomedcentral.com/1471-2105/8/S6/S5] 10.1186/1471-2105-8-S6-S5

    Article  Google Scholar 

  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): 6286-6291. 10.1073/pnas.0913357107.http://www.pnas.org/content/107/14/6286.abstract, 10.1073/pnas.0913357107

    Article  CAS  Google Scholar 

  23. Meinshausen N, Bühlmann P: High dimensional graphs and variable selection with the Lasso. Ann Stat. 2006, 34: 1436-1462. 10.1214/009053606000000281.

    Article  Google Scholar 

  24. Efron B, Hastie T, Johnstone I, Tibshirani R: Least angle regression. Ann. Stat. 2004, 32 (2): 407-499. 10.1214/009053604000000067.

    Article  Google Scholar 

  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, 33-40.

    Google Scholar 

  26. Meinshausen N, Bühlmann P: Stability selection. J R Stat Soc Ser B. 2010, 72 (4): 417-473. 10.1111/j.1467-9868.2010.00740.x.

    Article  Google Scholar 

  27. Tibshirani R: Regression shrinkage and selection via the lasso. J R Stat Soc Ser B. 1996, 58: 267-288.

    Google Scholar 

  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): 796-804. 10.1038/nmeth.2016.

    Article  CAS  Google Scholar 

  29. Breiman L: Random forests. Mach Learn. 2001, 45: 5-32. 10.1023/A:1010933404324.

    Article  Google Scholar 

  30. Weisberg S: Applied linear regression. 1981, New-York, Wiley

    Google Scholar 

  31. Hastie T, Tibshirani R, Friedman J: The elements of statistical learning: data mining, inference, and prediction. 2001

    Book  Google Scholar 

  32. Mairal J, Bach F, Ponce J, Sapiro G: Online Learning for Matrix Factorization and Sparse Coding. J Mach Learn Res. 2010, 11: 19-60. [http://jmlr.csail.mit.edu/papers/v11/mairal10a.html]

    Google Scholar 

  33. Schaffter T, Marbach D, Floreano D: GeneNetWeaver: in silico benchmark generation and performance profiling of network inference methods. Bioinformatics. 2011, 27 (16): 2263-2270. 10.1093/bioinformatics/btr373. [http://bioinformatics.oxfordjournals.org/content/27/16/2263.abstract] 10.1093/bioinformatics/btr373

    Article  CAS  Google Scholar 

  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): 229-239. 10.1089/cmb.2008.09TT. [http://online.liebertpub.com/doi/abs/10.1089/cmb.2008.09TT] 10.1089/cmb.2008.09TT

    Article  CAS  Google Scholar 

  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—D870-10.1093/nar/gkm815.

    Google Scholar 

  36. Gama-Castro S, Salgado H, Peralta-Gil M, Santos-Zavaleta A, Muñiz-Rascado L, Solano-Lira H, Jimenez-Jacinto V, Weiss V, García-Sotelo JS, López-Fuentes A, Porrón-Sotelo L, Alquicira-Hernández S, Medina-Rivera A, Martínez-Flores I, Alquicira-Hernández K, Martínez-Adame R, Bonavides-Martínez C, Miranda-Ríos J, Huerta AM, Mendoza-Vargas A, Collado-Torres L, Taboada B, Vega-Alvarado L, Olvera M, Olvera L, Grande R, Morett E, Collado-Vides J: RegulonDB version 7.0: transcriptional regulation of Escherichia coli K-12 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]

    Google Scholar 

  37. Küffner R, Petri T, Tavakkolkhah P, Windhager L, Zimmer R: Inferring gene regulatory networks by ANOVA. Bioinformatics. 2012, 28 (10): 1376-1382. 10.1093/bioinformatics/bts143.

    Article  Google Scholar 

  38. Mordelet F, Vert JP: SIRENE: Supervised inference of regulatory networks. Bioinformatics. 2008, 24 (16): i76—i82-10.1093/bioinformatics/btn273.

    Article  Google Scholar 

Download references

Acknowledgements

JPV was supported by the French National Reseach Agency (ANR-09-BLAN-0051-04 and ANR-11-BSV2-013-02) and the European Research Council (SMAC-ERC-280032).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Jean-Philippe Vert.

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

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Haury, AC., Mordelet, F., Vera-Licona, P. et al. TIGRESS: Trustful Inference of Gene REgulation using Stability Selection. BMC Syst Biol 6, 145 (2012). https://doi.org/10.1186/1752-0509-6-145

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1752-0509-6-145

Keywords