Bagging Statistical Network Inference from Large-Scale Gene Expression Data

Modern biology and medicine aim at hunting molecular and cellular causes of biological functions and diseases. Gene regulatory networks (GRN) inferred from gene expression data are considered an important aid for this research by providing a map of molecular interactions. Hence, GRNs have the potential enabling and enhancing basic as well as applied research in the life sciences. In this paper, we introduce a new method called BC3NET for inferring causal gene regulatory networks from large-scale gene expression data. BC3NET is an ensemble method that is based on bagging the C3NET algorithm, which means it corresponds to a Bayesian approach with noninformative priors. In this study we demonstrate for a variety of simulated and biological gene expression data from S. cerevisiae that BC3NET is an important enhancement over other inference methods that is capable of capturing biochemical interactions from transcription regulation and protein-protein interaction sensibly. An implementation of BC3NET is freely available as an R package from the CRAN repository.


Introduction
Gene networks represent the blueprint of the causal interplay between genes and their products on all molecular levels [1][2][3][4][5][6]. Gene regulatory networks (GRN) inferred from large-scale gene expression data aim to represent signals from these different levels of the gene network. The inference, analysis and interpretation of a GRN is a daunting task due to the fact that the concentrations of mRNAs provide only indirect information about interactions occurring between genes and their gene products (e.g., protein interactions). The reason for this is that DNA microarrays measure only the concentration of mRNAs rather than the binding, e.g., between proteins or between a transcription factor and the DNA. Despite the increased community effort in recent years [7,8] and a considerable number of suggested inference methods [9][10][11][12][13][14][15][16][17][18][19] there is an urgent need to further advance our current methods to provide reliable and efficient procedures for analyzing the increasing amount of data from biological, biomedical and clinical studies [20][21][22]. For this reason, this field is currently vastly expanding. A detailed review for many of the most widely used methods can be found in [15,18,[23][24][25][26].
A major problem for the inference of regulatory networks are the intricate characteristics of gene expression data. These data are high-dimensional, in the order of the genome size of the studied organism, and nonlinear due to the intertwined connection of the underlying complex regulatory machinery including the multilevel regulation structures (DNA, mRNA, protein, protein complexes, pathways) and turnover rates of the measured mRNAs, products and proteins. Further, gene expression data for network inference are large-scale, although, the ''Large p Small n'' [27] problem holds, because the number of explanatory variables (p genes) exceeds the number of observations (n microarray samples). In addition, technical noise and outliers can make it difficult to gain access to the true biological signal of the expression measurement itself.
The main contribution of this paper is to introduce a new network inference method for gene expression data. The principle idea of our method is based on bootstrap aggregation [28,29], briefly called bagging, in order to create an ensemble version of the network inference method C3NET [9]. For this reason we call our new method bagging C3NET (BC3NET). The underlying procedure of BC3NET is to generate an ensemble of bootstrap datasets from which an ensemble of networks is inferred by using C3NET. Then the obtained inferred networks are aggregated resulting in the final network. For the last step we employ statistical hypotheses tests removing the need to select a threshold parameter manually. Instead, a significance level with a clear statistical interpretation needs to be selected. This is in contrast with other studies, e.g., [30].
Given the challenging properties of gene expression data, briefly outlined above, BC3NET is designed to target these in the following way. First, BC3NET is based on statistical estimators for mutual information values capable of capturing nonlinearities in the data. Second, in order to cope with noise and outliers in expression data, we employ bagging because it has the desirable ability to reduce the variance of estimates [28]. Computationally, this introduces an additional burden, and a necessary prerequisite for any method to be used in combination with bagging is its tractability to be applicable to a bootstrap ensemble. C3NET is computationally efficient to enable this, even for high-dimensional massive data.
There are a few network inference methods that are similar to BC3NET. The method GENIE3, which was best performer in the DREAM4 In Silico Multifactorial challenge [31], employs also an ensemble approach, however, in combination with regression trees, e.g., in form of Random Forests [32]. In [30] a bootstrap approach has been used in combination with Bayesian networks to estimate confidence levels for features. However, we want to emphasize that, in contrast to BC3NET, both methods [30,31] do not provide a statistical procedure for determining an optimal confidence threshold parameter. Finally, we note that also for ARACNe a bootstrap version has been introduced [33], which has so far been used for inferring subnetworks around selected transcription factors [34,35].

Methods
The BC3NET approach for GRN inference In general, mutual information based gene regulatory network inference methods consists of three major steps. In the first step, a mutual information matrix is obtained based on mutual information estimates for all possible gene pairs in a gene expression data set. In the second step, a hypothesis test is performed for each mutual information value estimate. Finally, in the third step, a gene regulatory network is inferred from the significant mutual information values, according to a method specific procedure.
The basic idea of BC3NET is to generate from one dataset D(s), consisting of s samples, an ensemble of B independent bootstrap datasets fD b k g B k~1 by sampling from D(s) with replacement by using a non-parametric bootstrap [36] with B~1000. Then, for each generated data set D b k in the ensemble, a network G b k is inferred by using C3NET [9]. From the ensemble of networks fG b k g B k~1 we construct one weighted network which is used to determine the statistical significance of the connection between gene pairs. This results in the final binary, undirected network G. Fig. 1 shows a schematic visualization of this procedure. A base component of BC3NET is the inference method C3NET introduced in [9], which we present in the following in a modified form to obtain a more efficient implementation. Briefly, C3NET consists of three main steps. First, mutual information values among all gene pairs are estimated. Second, an extremal selection strategy is applied allowing each of the p genes in a given dataset to contribute at most one edge to the inferred network. That means we need to test only p different hypotheses and not p(p{1)=2. This potential edge corresponds to the hypothesis test that needs to be conducted for each of the p genes. Third, a multiple testing procedure is applied to control the type one error. In the above described context, this results in a network G b k . In order to test the statistical significance of the connection between gene pairs BC3NET utilizes the edge weights of the aggregated network G b w as test statistics. The edge weights of G b w are componentwise defined by Here I() is the indicator function which is 1 if its argument is G b k (i,j)~1 and 0 otherwise. This expression corresponds to the number of networks in fG b k g B k~1 which have an edge between gene i and j. For brevity, we write in the following n ij~G b w (i,j). From Eqn. 2 follows that n ij assumes integer values in f0, . . . ,Bg. Based on the test statistic n ij , we formulate the following null hypothesis which we test for each gene pair (i,j). with an edge between gene i and j is less than n 0 (a).
Here the cut-off value n 0 depends on the significance level a. Due to the independence of the bootstrap datasets we assume the null distribution of n ij to follow a binomially distributed Bin (B,p c ), whereas B corresponds to the size of the bootstrap ensemble and p c is the probability that two genes are connected by chance. The parameter p c relates to a population of networks, estimated from randomized data by using BC3NET, and corresponds to the fraction of randomly inferred edges in the bootstrap population (E½E b (B,D(s))) divided by the total number of possible edges in this population (E t (B)) that means The maximal number of gene pairs that can be formed from p genes in B bootstrap datasets is given by Figure 1. BC3NET algorithm: The gene regulatory network G is inferred from a bootstrap ensemble generated from a single gene expression dataset D. For each generated dataset in the ensemble, D b k , a network, G b k , is inferred using C3NET. From fG b k g B k~1 an aggregated network G b w is obtained whose edges are used as test statistics to obtain the final network G.
This value is independent of the sample size. E½E b (B,D(s)) corresponds to the expectation value of the number of randomly inferred edges for a population of an ensemble of bootstrap datasets of size B. Because E b (B,D(s)) is a random variable it is necessary to average over all possible bootstrap datasets of size B with sample size s. On a theoretical note we remark that these bootstrap datasets constitute a population that specifies a probability mass function (pmf) for which the expectation of E b (B,D(s)) needs to be evaluated. Due to the fact that this pmf is unknown the value of E½E b (B,D(s)) needs to be estimated. In order to estimate E½E b (B,D(s)) we randomize the data to estimate the number of edges randomly inferred in an bootstrap Using E b &E½E b (B,D(s)) as plug-in estimator for Eqn. 3 we obtain an estimate for p c . This allows us to calculate a p-value for each gene pair (i,j) and a given test statistic n ij , given by Eqn. 2, from the null distribution of n ij by Here p(i,j) is the probability to observe n ij or more edges by chance in a bootstrap ensemble of size B and sample size s.
Because we need to test p(1{p)=2 hypotheses simultaneously (one for each gene pair) we need to apply a multiple testing correction (MTC) [37,38]. For our analysis we are using a Bonferroni procedure for a strong control of the family-wise error rate (FWER). Typically, procedures controlling the FWER are more conservative than procedures controlling, e.g., the false discovery rate (FDR) by making only mild assumptions about the underlying data [39,40]. Based on these hypotheses tests the final network G is componentwise defined by That means if the connection between a gene pair is statistically significant they are connected by and edge, otherwise there is no connection.

Null-distribution of mutual information values
In order to determine the statistical significance of the mutual information values between genes we test for each pair of genes the following null hypothesis.
H I 0 : The mutual information between gene i and j is zero. Because we are using a nonparametric test we need to obtain the corresponding null distribution for H I 0 from a randomization of the data. Principally, there are several ways to perform such a randomization which conform with the formulated null hypothesis. For this reason, we perform 3 different randomizations and compare the obtained results with respect to the performance of the inference method to select the most appropriate one. Two randomization schemes (RM1 and RM2) permute the expression profiles for each gene pair separately. RM1 permutes only the sample labels and RM2 permutes the sample and the gene labels. In contrast, the randomization scheme RM3 permutes the sample and gene labels for all genes of the entire expression matrix at once.

Mutual Information Estimators
Due to the expected nonlinearities in the data we use mutual information estimators to assess the similarity between gene profiles instead of correlation coefficients. In a previous study, we found that for normalized microarray data the distribution among individual gene pairs can strongly deviate from a normal distribution [41]. This makes it challenging to judge by theoretical considerations only which statistical estimator is most appropriate for gene expression data because most estimators were designed assuming normal data. For this reason we compare eight different estimators and investigate their influence on the performance of C3NET.
Mutual Information is frequently estimated from the marginal and joint entropy H of two discretized random variables X and Y [42], In our study, we use four MI estimators based on continuous data and four MI estimators based on discretized data. The MI estimators for discretized data are the empirical estimator [42], Miller-Madow [42], shrinkage [43] and the Schürmann-Grassberger [44] mutual information estimator. For the emipirical estimator, the entropy H emp is estimated from the observed cell frequencies for each bin k of a random variable discretized into p bins, i.e., With an increasing number of bins, the empirical estimator underestimates the true entropy H due to undersampling of the cell frequencies n k . The different estimators attemp to adjust the undersampling bias by a constant factor [42], estimate cell frequencies by a shrinkage function between two models [43] or add a pseudo count from a probability distribution to the cell frequencies [44]. Mutual information can also be estimated from continuous random variables. The B-spline estimator considers the bias induced by the discretization for values falling close to the boundaries of a bin. For each bin, weights are estimated for the corresponding values from overlapping polynomial B-spline functions [45]. Hence, this method allows to map values to more than one bin.
For normal data, there is an analytical correspondence between a correlation coefficient and the mutual information [25], In this equation, the coefficient r could be the Pearson correlation coefficient r, Spearman rank correlation coefficient r or the Kendal rank correlation coefficient t.

Yeast gene expression data
We use the S. cerevisiae Affymetrix ygs98 RMA normalized gene expression compendium available from the Many Microbe Microarrays Database M3D [46]. The yeast compendium dataset comprises 9335 probesets and 904 samples from experimental and observational data from anaerobic and aerobic growth conditions, gene knockout and drug perturbation experiments. We map the yeast affymetrix probeset IDs to gene symbols using the annotation of the ygs98.db Bioconductor package. Multiple probesets for the same gene are summarized by the median expression value. The resulting expression matrix comprises a total of 9163 features for 4837 gene symbols and 4326 probesets that cannot be assigned to a gene symbol.

Simulated gene expression data
We simulate a variety of different gene expression datasets for Erdös-Rényi networks [47] with an edge density of f0:003,0:006,0:008,0:010g. An Erdös-Rényi network is generated by starting with n unconnected vertices. Then, between each vertex pair an edge is included with a pre-selected probability. The generated networks contain 150 genes of which f60,22,19,10g genes are unconnected. For each network, simulated gene expression datasets were created for various sample sizes of f50,100,200,500,1000g by using Syntren [48] including biological noise. We generate also simulated gene expression datasets for different subnetworks from the E.coli transcriptional regulatory network obtained from RegulonDB. The giant connected component (GCC) of the transcriptional regulatory network of E.coli consists of 1192 genes. We sample seven connected subnetworks from the GCC of sizes f50,100,150,200,300,400,500g. Again, using Syntren we simulate 100 different expression datasets including biological noise with sample size n~50 for each of these seven networks.

Gene pair enrichment analysis (GPEA)
To test the enrichment of GO-terms in the inferred yeast BC3NET network we adopt a hypergeometric test (one-sided Fisher exact test) for edges (gene pairs) instead of genes in the following way. For p genes there is a total of N~p(p{1)=2 different gene pairs. If there are p GO genes for a given GO-term then the total number of gene pairs is m~p GO (p GO {1)=2. Suppose the inferred yeast BC3NET network contains n edges of which k are among genes from the given GO-term, then a p-value for the enrichment of this GO-term can be calculated from a hypergeometric distribution by Here the p-value estimates the probability to observe k or more edges between genes from the given GO-term. For all GO:0032991 (macromolecular complex) offspring terms from Cellular Component that correspond to protein complexes, the above null hypothesis reflects the expected connection in a protein complex which is a clique (fully connected). For all other GO categories that we test, e.g., from the category Biological Process, the above is a very conservative assumption.

Influence of the randomization and MTC
The influence of the randomization scheme on the performance of BC3NET is shown in Fig. 2. Here we use simulated data from a Erdös-Rényi network consisting of 150 genes, of which 60 are unconnected. The figure shows results for RM1-RM3 with and without MTC for five different sample sizes, shown in the legend of the figure. As one can see, all three randomization schemes with a Bonferroni correction perform similarly good. Also RM1-RM3 without MTC perform similarly, however, significantly worse indicating the importance to correct for multiple hypotheses testing. Due to the fact that RM3 is from a computational point of view more efficient than RM1 or RM2 we use this randomization scheme for our following investigations. Fig. 2 includes also the F-scores obtained from the randomization of the expression data itself (right-hand side) to obtain baseline values for a comparison with the results from RM1-RM3. This is interesting because, e.g., in contrast to the AU-ROC [49], the F-score for data containing only noise is not 0:5 as for the AU-ROC. From this perspective, one can see that even the results without MTC are significantly better than expected by chance.

Influence of the mutual information estimator
To study the influence of the statistical estimators of the mutual information values, we use simulated data for several different network topologies. Fig. 3 shows results for eight different estimators and different sample sizes for a Erdös-Rényi network with an edge density of e~0:006. The three continuous estimators, Pearson, Spearman and Kendall as well as B-spline, perform better for smaller sample sizes. For large sample sizes the empirical, Miller-Madow, shrinkage and Schürmann-Grassberger perform slightly better. We want to note that for different parameters of the Erdös-Rényi network and different network types we obtain qualitatively similar results (not shown). Considering the size of the studied networks we used for our analysis, which contain 150 genes, sample sizes up to 100 lead to a realistic ratio of n=p * v 0:66 which one can also find for real microarray data. Larger ratios are currently and the near future hard to achieve. For this reason, we assess the results for smaller sample sizes as more important, due to their increased relevance for practical applications. Based on these results we use for the following studies the B-spline estimator.

Comparative analysis of BC3NET
Computational complexity. In [9] the computational complexity of C3NET has been estimated as O(n 2 ), where n corresponds to the number of genes. For BC3NET this means that its computational complexity is O(B|n 2 ). Here B is the number of bootstraps. In order to provide a practical impression for the meaning of these numbers, we compare the computational complexity between the ARACNe bootstrap network approach, described in [33], and BC3NET. We performed an analysis for a gene expression data set with 5000 genes and 200 samples. The ARACNe algorithm needed 22 hours for a single run that means to analysis one bootstrap data set. This results in a total time of 2200 hours (100|22 hours) for 100 bootstraps, which are about *92 days. In contrast, the BC3NET algorithm completed this task in only 28 minutes for all 100 bootstraps.
Comparative analysis using simulated data. In order to gain insight into the quality of BC3NET we study it comparatively by contrasting its performance with GENIE3 and C3NET. In Fig. 4 we show results for three different Erdös-Rényi networks each with 150 genes, of which f22,19,10g genes are unconnected. The edge density of these networks is [f0:003,0:006,0:008g. We use these edge densities because regulatory networks are known to be sparsely connected [50]. The F-score distributions for all studied conditions are larger for BC3NET. We repeated the above simulations for subnetworks from the transcriptional regulatory network of E. coli and obtained qualitatively similar results. This demonstrates the robustness of the results with respect to different network types and network parameters.
To emphasize the actual gain in the number of true positive edges with respect to C3NET, on which BC3NET is based, we present in Fig. 5 the percentage of the increase of inferred true positive edges for various network sizes ranging from 50 to 500 genes for subnetworks from E. coli. For the results shown in the left figure, we use a fixed sample size of n~50 and for the right figure the sample size equals the number of genes, i.e., n~p. For a fixed sample size (n~50) the BC3NET networks show an increase of true positives edges w45%, with a more prominent increase for the larger networks. Quantitatively, this observation is confirmed by a linear regression which gives a none vanishing positive slope of 0:039 and an intercept of 44:24%. Both parameters are highly significant with p-values * v10 {16 . For the datasets with variable sample sizes (n~p) the percentage of inferred true positive edges remains constant with an increasing network size and is around 30%. We want to note that the results for n~p assess the asymptotic behavior of BC3NET because the number of samples n increases linearly with the number of genes p. That means, asymptotically, the gain of BC3NET over C3NET is expected to be *30%. On the other hand, for real data for which pwn holds, the expected gain is much larger, as one can see from the left figure, reaching 70%.

Analysis of the regulatory network of yeast
Using BC3NET, we infer a regulatory network for a large-scale gene expression dataset of Saccharomyces cerevisiae. Due to the fact that for Saccharomyces cerevisiae no gold standard reference network is available to assess the quality of the inferred GRN we evaluate the resulting network by using functional gene annotations and experimentally validated protein interactions.
The yeast network inferred by BC3NET is a connected network that contains 9,163 genes and 27,493 edges with an edge density of 0:00065. The degree distribution of this network follows a power-law distribution, *k {a , with an exponent of a~4:38. We tested a total of 2159 GO-terms from the category Biological Process, whereas each GO-term contains less than 1000 annotated genes. From these, 525 (24:31%) test significant using a Bonferroni procedure indicating an enrichment of gene pairs for the corresponding GO-terms. The strongest enrichment of gene pairs we find in our analysis are for ribosome biogenesis, ncRNA and rRNA processing, mitochondirial organization, metabolic and catabolic processes and cell cycle. See Table 1 for an overview of the top 30 results.
One of the most reliable to detect (biochemical) interaction types that can be experimentally tested and that correspond to causal interactions, are protein-protein interactions from protein complexes. The reason therefore is that protein-protein interactions establish a direct connection between the proteins by forming physical bonds. Therefore we study the extend of protein complexes, as defined in the GO database [51], that are present in the yeast BC3NET network. We perform GPEA for 377 GOterms, which correspond to different protein complexes. From these we identify 94 protein complex terms with significantly enriched gene-pairs. The top 30 GO-terms of protein complexes we find are listed in Table 2. Some of the largest protein complexes detected in the BC3NET network are ribonucleopro- Figure 2. Influence of different randomization schemes (RM1, RM1 and RM3) and the multiple hypothesis testing correction on the network inference performance, measured by the F-score. The legend shows the used sample sizes. Each randomization scheme is used with and without a Bonferroni correction. The boxplots labeled 'random' correspond to randomly permuted data to get an impression for random Fscores. doi:10.1371/journal.pone.0033624.g002 tein complexes (789 edges) including the cytosolic ribosome (315 edges) and mitochondrial ribosome (142 edges). Further protein complexes present in the yeast BC3NET network are the proteasome complex (77 edges), proton-transporting ATP synthase complex (13 edges) and DNA-directed RNA polymerase complex (16 edges).
Finally, we study experimentally evaluated protein-protein interactions extracted from the BioGrid database (release 3:1:77) [52] and compare them with our yeast BC3NET network. First, we find that the yeast PPI network from BioGrid and our yeast BC3NET network have 4,723 genes in common. Further, we find a total of 878 BioGrid interactions among 1,043 genes that are present in the yeast BC3NET network. These interactions are distributed over a total of 282 separate network components, each consisting of 2 or more genes. Among these, we find 11 network components with a significant component size, where the largest significant component includes 147 genes and the smallest significant component includes 9 genes. Significance was identified from gene-label randomized data generating a null distribution for the size of connected network components of the 1,043 genes. The resulting p-values were Bonferroni corrected. For each BioGRID component that is nested in the yeast BC3NET network, we conduct a GO enrichment analysis. From this analysis we use the GO-term with the highest enrichment value to annotate the individual network components, see Table 3.
One of the most extensively studied biological processes in Saccharomyces cerevisiae is the cell cycle. For cell cycle the GPEA gives a gene-pair enrichment p-value of 2:6e{55, see Table 1. In Fig. 6 we show the largest network component of the cell cycle inferred by BC3NET that includes 304 genes and 423 edges. From this network, 57 edges are confirmed in BioGrid (violet edges), 25 edges are from protein complex units (GO) (green edges) and 7 edges are present in both databases (orange edges).

Discussion
From the analysis of BC3NET for the gene expression data set from S. cerevisiae, we find in addition to a significant enrichment of over 500 GO-terms in the category Biological Process, the significance of 94 GO-terms in Cellular Component for protein complexes. The largest complexes we identified are the ribosome (p~1:9e{310) and proteasome protein complex (p~1:3{104). There are two main reasons why edges of these protein complexes are highly abundant in the yeast BC3NET network. First, the ribosome and proteasome protein complexes are well annotated because they have been extensively studied in yeast [53]. Second, the ribosome and proteasome protein complex are mainly regulated on the gene expression level and, where observed, having highly dependent gene expression patterns [53]. Therefore, it is plausible that GRN inference methods can also pick-up signals from physical interactions between protein subunits of protein complexes.
We want to note that we are not the first to recognize that gene expression data contain information about protein-protein interactions. For example, [53,54] provide evidence that proteins from the same complex show a significant coexpression of their corresponding genes. Also in [55] it is mentioned that inferred interactions from gene expression data 'may represent an expanded class of interactions' [55]. However, when it comes to the experimental assessment of the inferred networks, usually, only interactions related to the transcriptional regulation are studied, e.g., with ChIP-chip experiments [11,16]. To our knowledge we are the first to provide a large-scale analysis of an inferred GRN from gene expression data with respect to the presence of proteinprotein interactions.
BC3NET is an ensemble method that uses as base network inference algorithm C3NET [9,56]. As for other ensemble methods based on bagging, e.g., random forests, the interpretability and characteristics of the base method does usually not translate to the resulting ensemble method [28,32]. In our case this means that the inferred network can actually have more than p edges, despite the fact that networks inferred by C3NET can not. However, in our case this is a desirable property because it improves BC3NET leading ultimately to a richer connectivity structure of the inferred network. Specifically, our numerical results demonstrate that BC3NET gains in average more than 40% true positive edges compared to C3NET (see Fig. 5). Another more general advantage of an ensemble approach is that it is straight forward to use on a computer cluster because a parallelization is naturally given by the base inference methods. Given the increasing availability of computer clusters this appears to be a conceptual advantage over none ensemble methods, likely to gain even more importance in the future. In this paper we pursued a conservative approach by using a Bonferroni procedure for MTC to demonstrate that even in this setting our method is capable of inferring many significant interactions that can be confirmed biologically. However, there is certainly potential to use more adopted MTC procedures that are less conservative. For example, procedures controlling the false discovery rate (FDR) could be investigated [39,57].
Further, we want to note that despite the fact that the network inference method C3NET is no Bayesian method [58,59], BC3NET is. The reason for this is that it is known for the bootstrap distribution of a parameter to correspond approximately to the Bayesian posterior distribution for a noninformative prior, and the bagged estimate thereof is the approximate mean of the Bayesian posterior [60]. Hence, BC3NET can be considered as a Bayesian method with noninformative priors for the connectivity structure among the genes. Given the problem to define informative priors for a Bayesian approach in a genomics context, either because not enough reliable information about a specific organism is available or because it is difficult to select this information in an uncontroversial manner, a noninformative prior is in the current state of genomics research still a prevalent choice. From a theoretical point of view, a bootstrap implementation is false negative edges. However, as demonstrated by our numerical analysis, BC3NET is an important improvement toward the inference of causal gene regulatory networks.
Despite the fact that the presented inference method BC3NET was introduced by using gene expression data from DNA microarray experiments, it can also be used in connection with data from RNA-seq experiments. Given the rapidly increasing importance of this new technology we expect that within the next few years datasets with sufficient large sample size are available to infer GRN.