TIGRESS: Trustful Inference of Gene REgulation using Stability Selection

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 regulations has many potential far-reaching applications in biology and medicine, ranging from the in silico modelling and simulation of the gene regulatory network (GRN) to the identification of new potential drug targets.However, while many TF-TG regulations 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 [2,21,9,1,36,33,16,8,5,4]. 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 or mutual information [7,26,11] or take explicit advantage of perturbations in the experiments such as gene knock-downs [31].The difficulty to separate direct from indirect regulations has been addressed with the formalism of Bayesian networks [14,17,13], or by formulating the GRN inference problem as a feature selection problem [19].We refer to [27,24] 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 [28].This idea underlies the Bayesian network formalism [14], but is more directly addressed by GENIE3 [19], 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 [19,24].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) [10] combined with stability selection [3,29], for GRN inference.LARS is a computationally efficient procedure for multivariate feature selection, closely related to Lasso regression [34].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 [23].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.Overall this study confirms the potential of state-of-the-art feature selection techniques for GRN inference.

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 regulations of the form (t, g) for t ∈ T and g ∈ G.
We do not try to infer self-regulation, meaning that for each target gene g ∈ G we define the set of possible regulators as T g = T \{g} if g ∈ T 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 : E → R 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.

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 [7,26,11].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 post-process the predicted regulations and try to remove regulations likely to be indirect because they are already explained by other regulations [26].Another strategy is, given a target gene g ∈ G, 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 [28,19].More specifically, for each target gene g ∈ 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 ∈ T g : where X i represents the expression level of the i-th gene across different experiments (modelled as a random variable), X Tg = {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 then the score s g (t) should typically assess the probability that β t,g is non-zero [28].More general models are possible, for example [19] model f g with a random forest [6] 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 [10].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 [35,18], 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 [34,18], 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, [10]) and MATLAB (SPAM toolbox, [22]).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 [29] 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 [3,29].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/2 times 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 > 0 steps 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 ∈ G, t ∈ T g and l ∈ [1, L] the 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 divide the frequency by R to 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∈Tg F (g, t, l) = l, for any gene g and LARS step l.   , 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 [3,29] 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: 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 = 1 meaning 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 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 g ∈ G 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 α = 1 means that no weight randomization is performed on the different expression arrays, while α = 0 leads to maximal randomization.[29] advocate that a value between 0.2 and 0.8 is often a good choice.Regarding the choice of L, [29] 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.

Other GRN inference methods
We experimentally compare TIGRESS to several other GRN inference methods.We use the MATLAB implementations of CLR [11] and GENIE3 [19].We run ARACNE [26] using the R package minet.We keep default parameter values for each of these methods.Results borrowed from the DREAM5 challenge [23] were directly obtained by each participating team.

Performance evaluation
Given a gene expression data matrix, each GRN inference method outputs a ranked list of putative TF-TG regulations.Taking only the top K predictions in this list, we can compare them to known regulations to assess the number of true positives (T P , the number of known regulations in the top K predictions), false positives (F P , the number of predicted regulations in the top K which are not known regulations), false negatives (F N , the number of known interactions which are not in the top K predictions) and true negatives (T N , 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 (T P/(T P + F P )), recall (T P/(T P + F N )), and fall-out (F P/(F P + T N )).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 [23] for more details on this scoring scheme.Finally, we compute the so-called overall score for a GRN inference method by integrating the AUROC and AUPR P -values as follows: 3 Data We evaluate the performance of TIGRESS and other GRN inference methods on four benchmark datasets, each consisting of a compendium of gene expression data, a list of known TF, and a gold standard set of verified TF-TG regulations 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 four networks.
The first three benchmarks are taken from the DREAM5 challenge [23].Network 1 is a simulated dataset.Its topology and dynamics were modelled according to known GRN, and the expression data were simulated using the GeneNetWeaver software [32].We refer the interested reader to [25,24] 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 is provided for this dataset consisting in expression data for S. aureus.
Additionally, we run experiments on the E. coli dataset from [11], 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 [12] (version 4 build 6).It consists in 907 experiments and 4297 genes.We obtained the gold standard data from RegulonDB [15] (version 7.2, May 6th, 2011) that contains 3812 verified interactions among 1525 of the genes present in the microarrays experiments.
As a pre-processing step, we simply mean-center and scale to unit variance the expression levels of each gene within each compendium.

DREAM5 challenge results
In 2010 we participated to the DREAM5 Network Inference Challenge, an open competition to assess the performance of GRN methods [23].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 overall score (5), and an overall ranking was proposed based on the average of the overall scores over the three networks.We submitted predictions for all networks with a version of TIGRESS which we did not try to optimize, which we refer to 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 single-core 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.2 followed previous suggestions for the use of stability selection [29], 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 in average overall score.The winning method, both in silico and overall, was the GENIE3 method of [19].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 [20] ranked second overall, with particularly strong performance on 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 TF-TG regulation.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 TI-GRESS 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.2 and 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}.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.2 and 0.8, and any L less than 10 leads to an overall 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 [29] 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.3 and 0.5, suggesting that the adjustment to the number of bootstraps can mostly be made through the choice of L. Finally, we unsurprisingly observe that increasing the number R of resampling runs leads to better performances.On Figure 4, we show the overall 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.  3 summarizes these performances in terms of AU P R, AU ROC and related p-values as well as the overall score.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.

Comparison with other methods
TIGRESS outperforms all methods in terms of AUPR and all methods but GENIE 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   The performance of GENIE3, Naive TIGRESS and ANOVA were obtained during the DREAM5 competition.TIGRESS corresponds to the choice of parameters leading to the best performance (area score, α = 0.4, L = 2, R = 8, 000).We ran CLR and ARACNE using public implementations of these methods.

Method
TIGRESS is more reliable than GENIE in its first predictions but contains overall more errors when we go further in the list.

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).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 vitro networks.As a partial answer, Figure 6 shows the behavior of the overall score with respect to L for Networks 3 and 4. Interestingly, it seems that a larger L is preferable in this case, suggesting that one may have to adapt the parameters to the size of the network.Indeed, networks 3 and 4 contain respectively 4, 511 and 5, 950 nodes, making them much larger than the in silico network we tuned the parameters on.However, the improvement is not dramatic in absolute value.
On Figure 7 we compare Precision/Recall and ROC curves obtained with TIGRESS with several other algorithms on the E. coli network from [11].Table 5 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 analyse the type of errors it typically makes taking the E. coli dataset as example.We analyse 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 8, 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 px the proportion of spurious TG-TF couples with distance x.  Figure 9 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 distribution H(N S , px 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 distance-2 edges in a significantly higher than p2 whereas significantly less distance-> 4 edges 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 Table 6, we detail the three possible patterns observable in this situation.

Method
Figure 10 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.Grandparent/Grandchild G1 has a child that is a parent of G2.G1 is a grandparent of G2.

Discussion
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 [29] 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 feature t can be both expressed in a common formalism as E [φ(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 the in silico network, we observed that it achieves a similar performance to that of GENIE3 [19], the best performer at the DREAM5 challenge.However, TIGRESS does not do as good as this algorithm on in 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 [11] showed that many spuriously detected edges follow the same pattern: TIGRESS discovers edges when in reality the two nodes are siblings, 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 as switches, 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 that in vivo networks gold standards may not be complete.Therefore, the hypothesis that TIGRESS is actually correct when predicting these loops cannot be discarded.
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 on in vivo networks.For example, performances on the E. 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 in the global approach such as adding some supervision in the learning process as, e.g., investigated in [30].

Figure 1 :
Figure1: 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.

Figure 2
Figure 2 summarizes the overall 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.2 and 0.8, and any L less than 10 leads to an overall 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[29] 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.

Figure 2 :
Figure 2: Overall score for Network 1. From top to bottom, plots show the results for R = 1, 000, R = 4, 000 and R = 10, 000 for both the area (left) and the original (right) scoring settings, as a function of α and L.

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

Figure 5
Figure 5 depicts both the ROC and the Precision/Recall curves for several methods on Network 1. Table3summarizes these performances in terms of AU P R, AU ROC and related p-values as well as the overall score.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 outperforms all methods in terms of AUPR and all methods but GENIE 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

Figure 4 :
Figure 4: Overall score as a function of R. In both scoring settings, α and L were set to 0.4 and 2, respectively.

Figure 5 :
Figure 5: ROC (left) and Precision/Recall (right) curves for several methods on Network 1.

Figure 7 :
Figure 7: Precision/Recall (Left) and ROC (Right) curves for several methods on the E. coli dataset.

Figure 8 :
Figure 8: Exact distribution of the shortest path between spuriously predicted TF-TG couples.
discovered edges (by order of discovery)

Figure 9 :
Figure 9: 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.

Figure 10 :
Figure 10: Distribution of distance 2 errors.95% error bars were computed using the quantiles of a hypergeometric distribution.

Table 1 :
Datasets used in our experiments.

Table 2 :
AUPR, AUROC and minus the logarithm of related p-values for all DREAM5 Networks and the three best teams.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.

Table 3 :
Comparison of different methods on Network 1 of the DREAM5 challenge.

Table 4 :
Table 4 shows the values of AUPR, AUROC, related p-values and overall score for DREAM5 networks 3 and 4 reached by TIGRESS.TIGRESS performance on DREAM5 Networks 3 and 4.

Table 5 :
TIGRESS compared to state-of-the-art methods on the E. coli Network.

Table 6 :
Distance-2 patterns between two nodes G1 and G2 in a directed graph.