 Methodology article
 Open Access
 Published:
Fast Bayesian inference for gene regulatory networks using ScanBMA
BMC Systems Biologyvolume 8, Article number: 47 (2014)
Abstract
Background
Genomewide timeseries data provide a rich set of information for discovering gene regulatory relationships. As genomewide data for mammalian systems are being generated, it is critical to develop network inference methods that can handle tens of thousands of genes efficiently, provide a systematic framework for the integration of multiple data sources, and yield robust, accurate and compact genetogene relationships.
Results
We developed and applied ScanBMA, a Bayesian inference method that incorporates external information to improve the accuracy of the inferred network. In particular, we developed a new strategy to efficiently search the model space, applied data transformations to reduce the effect of spurious relationships, and adopted the gprior to guide the search for candidate regulators. Our method is highly computationally efficient, thus addressing the scalability issue with network inference. The method is implemented as the ScanBMA function in the networkBMA Bioconductor software package.
Conclusions
We compared ScanBMA to other popular methods using time series yeast data as well as timeseries simulated data from the DREAM competition. We found that ScanBMA produced more compact networks with a greater proportion of true positives than the competing methods. Specifically, ScanBMA generally produced more favorable areas under the ReceiverOperating Characteristic and PrecisionRecall curves than other regressionbased methods and mutualinformation based methods. In addition, ScanBMA is competitive with other network inference methods in terms of running time.
Background
Identifying gene regulatory networks is an important problem in biology. There have recently been many advances in this area in terms of tools for collecting and analyzing largescale genomics data. Many of these datasets, from microarrays and next generation sequencing, quantify the expression levels of all genes in a given genome. Genomewide timeseries data, in principle, allow reverse engineering of the gene regulatory relationships by studying the temporal patterns of regulators and target genes. However, this can be a difficult problem due to the large number of genes (i.e. variables) being measured, which typically far exceeds the number of observations. Also, the number of actual regulators for a particular gene is only a small fraction of the number of possible regulators.
A popular method for inferring gene regulatory networks from time series data uses Dynamic Bayesian Networks (DBN)[1–5]. DBN methods estimate a probabilistic graphical model, given the timeseries data. DBN methods work well, but the network size that they can handle in practice is limited because of their computational cost.
Ordinary differential equations (ODE) are alternative methods for constructing networks[6, 7]. These methods are deterministic rather than statistical, although ODE methods can be combined with statistical methods. DBN on local networks within a larger ODE model inference method have been used, for example[8].
Another class of methods is based on regression models in which parent nodes (regulators) are inferred for each target gene. Vector autoregressive models have been proposed for inferring causal links between genes. Often this takes the form of a model selection problem, and methods such as the Least Absolute Shrinkage and Selection Operator (LASSO)[9, 10], elastic net[11, 12], and Bayesian model averaging (BMA)[13, 14] have been used[15–19]. Morrissey, et al.[20] implemented a Markov Chain Monte Carlo (MCMC) sampler for a fully Bayesian formulation of the autoregressive model.
Mutual information methods have been used extensively on genetics data[21–24], but usually for steadystate rather than timeseries data. These methods are typically nondirectional. Recently, mutual information methods have been extended to analyze timeseries data and produce directed networks[25, 26]. Mutual information methods have the advantage of being able to identify nonlinear relationships.
Our contributions
We present a new approach using Bayesian Model Averaging (BMA) for variable selection from timeseries gene expression data. Our new method offers the following advances over our previous work[18, 19]:

We develop a new algorithm called ScanBMA that searches the model space more efficiently and thoroughly than previous algorithms. It is much faster than previous implementations of BMA for a large number of predictors, resulting in running time comparable to that of LASSO. It allows inference for networks of thousands of genes to be completed in minutes on a standard laptop computer.

We transform the timeseries data to reduce spurious correlations. Specifically, we remove the effect of a gene on itself by subtracting the mean expression level for each gene at each time point and then using the residuals from a regression of its expression at the current time point on its expression at the previous time point.

We use Zellner’s gprior[27] for the regression parameters and show that using the gprior to compute posterior probabilities outperforms our previous effort using the Bayesian Information Criterion (BIC).

We address the scalability of network inference methods. Our new implementation produces running times of minutes compared to hours or even days for some competing methods, thus offering substantial improvements.
We also carried out extensive empirical studies of our new method. Specifically, we compared our new method, ScanBMA, to other network construction methods in the literature, namely LASSO, the mutual information methods MRNET (Maximum Relevance/Minimum Redundancy)[24], CLR (Context Likelihood or Relatedness)[23] and ARACNE (Algorithm for the Reconstruction of Accurate Cellular Networks)[22], and also Dynamic Bayesian Networks (DBN) when the latter were computationally feasible.
We benchmarked the performance of our approach, ScanBMA, using two datasets. The first dataset measures the gene expression levels over time of 97 yeast segregants perturbed with the drug rapamycin. The second dataset consists of simulated timeseries data from the DREAM4 (Dialogue for Reverse Engineering Assessment and Methods) challenge. For the yeast dataset, we found that our method outperformed competitors and previous analyses in recovering regulatory relationships previously reported in the literature. For the DREAM4 data, for which no prior information was available, our method performed comparably to other methods, while producing more compact networks. Finally, the ScanBMA algorithm presents a substantial improvement in running time over previous implementations of BMA. The method is implemented as the ScanBMA function in the networkBMA Bioconductor software package.
Results and discussion
Method outline
In ScanBMA, network inference is formulated as a series of variable selection problems in which parent nodes (regulators) are inferred for each target gene. The BMA framework accounts for model uncertainty in variable selection by averaging over the posterior distributions from multiple models, weighted by their posterior model probabilities[13, 14]. A challenge of BMA is to efficiently select a set of models to be averaged over. ScanBMA uses a greedy approach to explore the model space and uses the Occam’s window principle[28] to eliminate unlikely models.
We previously developed a supervised framework to integrate external data sources, including coexpression, genomewide binding, sequence polymorphism, physical interaction, genetic interaction, and literature curation data[18]. Using a training set consisting of approximately 500 known regulatory relationships in the literature, we computed prior probabilities of regulatory relationships across all candidate genes and regulators. These prior probabilities were then used to compute the posterior probabilities of of candidate regression models. We used Zellner’s gprior[27] to specify the prior for the model parameters in ScanBMA. We developed an expectationmaximization (EM) algorithm to estimate the prior variance parameter g.
Before the regression step, we apply a univariate measure (such as Rsquared or BIC) to rank candidate regulators for each target gene using these prior probabilities of regulatory relationships. The parameter nvar controls the number of top regulators used in the regression step of each target gene. We have performed empirical studies to study the effect of and estimate the optimal nvar.
Assessment
A number of metrics have been used to evaluate the quality of inferred networks. We focus on a few that compare the inferred network with a goldstandard network of true edges. One measure that we use is the precision of the inferred network, equal to the number of true positives divided by the total number of edges in the inferred network. Precision is important to researchers because an experiment to verify relationships identified in exploratory work can be expensive. Thus, the more confident we can be when identifying relationships, the better. In light of this, we also look at the area under the precisionrecall curve (AUPRC). This gives a more comprehensive view of network quality and does not require that a threshold be chosen for the posterior probability of an edge or for the number of edges included. We also look at the area under the ROC curve (AUROC), which is widely used to assess the quality of networks.
Due to incomplete knowledge in real data, we use a partial assessment based on the YEASTRACT database[29]. This is a literaturecurated repository of regulatory relationships between known transcription factors and target genes in yeast, based on more than 1,300 literature references. Due to incomplete knowledge in yeast biology, the lack of an edge in YEASTRACT is not hard evidence of the absence of a relationship between two genes, although it is used as such in our evaluation. In contrast, the true underlying networks from the simulated DREAM4 data are known, so that an absence of a true edge is a false negative and a presence of a nonexisting edge is a false positive.
Results: yeast time series data
Table1 summarizes the assessment results for different methods applied to the yeast dataset. Our new ScanBMA method with nvar = 20 had a precision of 0.39, much higher than any other method, including our previous method, iBMA[19]. The area under the ROC curve (AUROC) for ScanBMA was also much better than for the competing methods. Note that for random guessing, the expected AUROC is 0.5 and the expected AUPRC is 0.038. For these data, the mutual information methods CLR and MRNET produced very large networks — nearly half of the number of possible edges, which is not concordant with what is known about regulatory networks of this type.
Table2 compares the preferred version of ScanBMA with four other versions, each of which lacks one of the components of our final method. This table shows that each component contributes to the accuracy of our final network. As shown by[19], the incorporation of the informative prior yields the largest increase in performance, but the other components also contribute. The use of the gprior reduces the number of false positives, while the data transformation substantially reduces the size of the inferred network.
The precisionrecall curves for the different methods on the yeast data are shown in Figure1. This figure shows the quality of the ScanBMA network even beyond the 95% posterior inclusion probability cutoff. The precision stays high through a large range of recall, whereas for the other methods it quickly drops to the level of random guessing. The figure also shows that nvar = 20 performs better than nvar = 3556.
When analyzing gene networks where the number of nodes is in the thousands, computation time can be an important consideration. We compared ScanBMA with the other methods on the yeast data by running each method on 20 target genes under controlled conditions to find the average cpu time per gene.
Table3 shows that ScanBMA with nvar = 3556, i.e. considering all other genes whose expression varied as potential regulators, is within a factor of 3 of LASSO, the fastest of the other methods. Some of the mutual information methods, on the other hand, are much slower, with MRNET taking about 50 times longer than ScanBMA. Table3 also shows that ScanBMA produces a substantial improvement in computational efficiency over our previous method iBMA[19], especially when nvar is large. Dynamic Bayesian network methods were not included in the comparison because they analyze the entire network at once and do not scale well enough to be run feasibly on large network datasets such as the yeast data.
Results: simulated DREAM4 data
Table4 summarizes the results of the competing methods for the DREAM4 10gene networks. For the DREAM4 networks, we were able to add the Dynamic Bayesian Networks method, as implemented in the ebdbnet R package[30], to the comparison because the networks are small. ScanBMA again performed best among the methods, particularly in terms of the areas under the ROC and precisionrecall curves, even though no external information was available. However, the extent of ScanBMA’s superiority to other methods, notably LASSO, was smaller in this case than for the yeast data, reflecting the lack of external information.
The precisionrecall curves from the various methods are shown for the first of the five 10gene networks in Figure2. Figure3 shows the actual networks returned by each method in comparison with the true network. The ScanBMA network resembles the true network fairly well. In particular, the compactness of the ScanBMA network is apparent, particularly when compared with LASSO and MRNET. The small number of false positives may be useful in focusing the attention of the biologist on edges of high interest when searching for new regulatory relationships.
The results of the methods for the DREAM4 100gene networks are summarized in Table5. For these networks, the mutual information methods MRNET and CLR performed best in terms of the area under curve measures. ScanBMA was not quite as good by these measures, but its precision was much higher than that of any other method. Figure4 illustrates the precisionrecall curves for the various methods, showing increased precision across a broader range for ScanBMA.
Conclusions
We have presented a Bayesian Model Averaging method for inferring gene regulatory networks from time series data. It incorporates external information in a principled way via the prior edge probabilities, transforms the data to reduce spurious correlations, and uses Zellner’s gprior for model parameters, with g estimated from the data. We have introduced a new algorithm, ScanBMA, to search the model space efficiently. Our method infers compact networks with higher precision than the competing methods we have assessed, important features for further analysis in searching for new regulatory relationships.
We found that our method outperformed previous methods as well as LASSO and mutual information methods on yeast timeseries data. In addition, our method performed comparably to competing methods, including Dynamic Bayesian Networks, on simulated data from the DREAM4 challenge, even in the absence of prior information. The networks from ScanBMA are also similar in size to the target networks.
Methods
Bayesian model averaging
In regressionbased methods for network inference, we infer regulators (parent nodes) for each target gene. Hence, network inference can be formulated as a series of variable selection problems. We use the following multiple linear regression model for network inference from timeseries data:
where X _{ i,t,s } is the expression level for gene i at time t for strain or replicate s, H is the set of potential regulators, and${\epsilon}_{i,t,s}\stackrel{\text{iid}}{\sim}N(0,{\sigma}_{\epsilon}^{2})$ is the error term, for i = 1,…,n, t = 2,…,T and s = 1,…,S. We are particularly interested in whether β _{ hi } ≠ 0, indicating a regulatory relationship from gene h to gene i.
One difficulty is that for gene network time series data there are typically far more potential regulators than observations. To address this problem, as well as to take into account uncertainty in model selection, we use BMA to obtain the posterior probability that each regulator is in the model[13, 14]. BMA takes model uncertainty into account by averaging over the posterior distributions of a quantity of interest based on multiple models, weighted by their posterior model probabilities. Thus the posterior probability that β _{ hi } ≠ 0, also called the posterior inclusion probability of β _{ hi }, is
where$\text{Pr}({M}_{k}D)\propto \int p(D{\theta}_{k},{M}_{k})p({\theta}_{k}{M}_{k})d{\theta}_{k}$ is the posterior probability of model M _{ k } given data D, p(Dθ _{ k },M _{ k }) is the likelihood of the parameter vector θ _{ k } of model M _{ k }, the prior probability of model M _{ k } is$\text{Pr}({M}_{k})\propto {\prod}_{h\in {M}_{k}}{\pi}_{\mathit{\text{hi}}}{\prod}_{h\notin {M}_{k}}(1{\pi}_{\mathit{\text{hi}}})$, π _{ hi } is the prior probability that gene h regulates gene i, and the models considered are denoted by M _{1},…,M _{ K }.
An additional issue is that the number of possible models is too large to enumerate in a reasonable amount of time. MCMC approaches have been applied where the number of potential regulators is large[31], but they are expensive computationally. Iterative BMA has also been applied to large genetics datasets, but the number of regulators was restricted prior to analysis[18, 32].
Data transformations
One concern when identifying regulatory relationships among genes is that there is a great deal of variation in gene expression levels that does not come from these interactions. For example, in the yeast data, many of the genes experience a sharp change in expression level over the first few time points caused by the application of the drug. This common trajectory is not important for inference and can produce many large correlations that do not correspond to actual interactions. In addition, we have found that removing the effect of a gene on itself can improve inference. By doing this we are removing excess variation in order to gain accuracy in inferring relationships.
In light of these observations, we transform the data in two steps to remove these extra sources of variation. First, we use timeadjusted data by subtracting the mean expression level for each gene i at timepoint t across strains:
where${\stackrel{\u0304}{X}}_{i,t,\xb7}={S}^{1}{\sum}_{s=1}^{S}{X}_{i,t,s}$. This removes the overall effect of the drug.
Second, we take the residuals from regressing the gene on itself at the previous timepoint:
where α _{ i } is the regression coefficient in the simple linear regression model of the expression level of gene i at time t on its expression level at time t  1:
where${\delta}_{i,t,s}\stackrel{\text{iid}}{\sim}N(0,{\sigma}_{i}^{2})$. Then${\widehat{\alpha}}_{i}$ is the least squares estimate of α _{ i } based on all S(T  1) relevant observations for gene i.
Specifying prior distributions
An important feature of our method is a way to incorporate external information. The Bayesian approach requires the specification of prior distributions and so includes this directly.
BMA requires two types of prior information: the prior edge probability, π _{ hi }, that potential regulator h regulates target gene i, for each h and i, and the prior distribution of the parameter vector for each model considered.
For the prior edge probabilities π _{ hi }, we considered two different prior distributions. The first is based on the empirical finding of Guelzim et al.[33] that each target gene is regulated by a small number of transcription factors, estimated as 2.76 per gene[19]. We implement this by setting π _{ hi } = 2.76/6000 for all h and i. This prior distribution does not incorporate any genespecific information. We call it the Guelzim prior.
The second prior edge distribution we consider is that of Lo et al.[19], which uses external information in the form of expression data, bindingsite data, curated literature and other sources to come up with a prior probability for each individual possible regulatory relationship. This leads to values of π _{ hi } specific to each regulatortarget pair. We refer to this as the informative prior. Integration of multiple information sources has been shown to be beneficial in network construction[34, 35].
As a prior for the model parameters, corresponding to the strengths of the relationships, we use Zellner’s gprior[27], as in[36]. The prior distribution of the parameter vector$({\beta}_{0,i}^{[k]},{\mathbf{\beta}}^{[k]},{\sigma}^{2[k]})$ of model M _{ k } is
where X _{ k } is the design matrix for model M _{ k } and g > 0 controls the prior variance of the regression parameters. This prior yields an analytic form for the posterior model probabilities P r(M _{ k }D), namely
where M _{0} is the null model with no regulators, d _{ k } is the number of regulators present in model M _{ k }, and${R}_{k}^{2}$ is the R ^{2} value for M _{ k }. This is an alternative to using BIC to approximate the posterior model probabilities, as has been done previously[32, 37].
The parameter g controls the expected size of the regression parameters β _{ hi }, and is approximately equal to the prior variance of${\beta}_{\mathit{\text{hi}}}/\text{SE}({\widehat{\beta}}_{\mathit{\text{hi}}})$, where${\widehat{\beta}}_{\mathit{\text{hi}}}$ is the OLS estimator. This suggests that g should be at least 1, otherwise the individual β _{ hi }’s would be expected to be nearly indistinguishable from the noise even under ideal conditions. Also, the effective number of data points in the prior is n/g. Using g = n corresponds to a unit information prior and yields similar results to using the BIC approximation. Raftery[38] argued that the prior should be no more spread out than the unit information prior, suggesting a choice of g in the range 1 ≤ g ≤ n.
We estimate g from the data by maximum marginal likelihood, where the likelihood is summed over the model space. We do this using an EM algorithm[39, 40], where the “missing data” are the γ _{ hi }’s indicating whether regulator gene h is in the model for target gene i. First, we run BMA with a selected value of g, and from this we obtain the posterior probabilities of the models used. We then maximize
with respect to g. We use the new value of g in the next iteration of ScanBMA. We do this until convergence.
Searching the model space using ScanBMA
To find the models to be included in the BMA summation, we use the Occam’s window principle, according to which models that are substantially worse than the best model found (by a factor of C in posterior probability) are discarded[28]. We used C = 100, based on an extensive review of conventional standards of scientific evidence. To find the models in Occam’s window, we propose a new, fast algorithm called ScanBMA.
The main idea of ScanBMA is to keep an active set of models around which to search. All models which can be created by adding or removing a single predictor from those in the active set are then evaluated and, if they fall within Occam’s window of the current best model, are added to the active set to expand the search. The method is outlined in Algorithm 1.
Algorithm 1: ScanBMA
Because ScanBMA does not average over every model in the model space, the algorithm yields posterior inclusion probabilities of either 100% or 0% for many regulators. These extreme posterior inclusion probabilities are only approximations, and we can refine them. To estimate the posterior inclusion probability of a predictor X _{ h } with an approximate posterior probability of 100%, we calculate the ratio O _{h,i } of the posterior probability of the best model to that of the best model with the predictor X _{ h } removed. We then approximate the posterior probability of predictor h by O _{h,i }/(1 + O _{h,i }), which will be less than 100%.
Similarly, to estimate the posterior inclusion probability of a predictor X _{ h } with approximate posterior probability 0%, we compute O _{ h,i }, the ratio of the posterior probability of the best model to that of the best model with the predictor X _{ h } added. Then the approximate posterior inclusion probability of predictor h is (1 + O _{ h,i })^{1}, which will be greater than 0%.
This postprocessing step yields a unique ordering among the predictors, which is useful in evaluation.
A variant of ScanBMA that we found to greatly improve computational efficiency without degrading performance is to restrict the number, nvar, of potential regulators h considered for each target gene i to those with the highest prior probabilities, π _{ hi }. Guelzim et al.[33] found that the number of transcription factors per gene has approximately an exponential distribution with mean 2.76; they did not find any of the target genes that they considered to have more than 13 regulators. As a result, we consider nvar = 20, which leaves some margin over the Guelzim maximum. We compare this with considering all genes with observed variation in expression as potential regulators, which amounts to setting nvar = 3556.
Data
To validate our method, we applied it to timeseries data from a gene expression experiment on yeast as well as simulated timeseries datasets from the DREAM4 competition.
The yeast data come from an experiment on 97 strains of yeast crossed from two parent strains[18]. Each strain was subjected to a treatment of the macrolide drug rapamycin, chosen to cause changes in gene transcription across the genome. Gene expression levels were measured for approximately 6,000 genes every 10 minutes from 0 to 50 minutes after the administration of the drug, using Affymetrix microarrays. The data were then filtered to remove genes that did not show significant variation over the measurements, leaving 3,556 genes.
The yeast data can be represented as a threedimensional array consisting of 3,556 genes, 97 segregants and 6 time points. There are no replicates in these data. However, the segregant axis and the time axis capture the genetic and temporal variations. These yeast data are publicly available from ArrayExpress with accession number EMTAB412.
The simulated DREAM4 In Silico Network Challenge[41–43] provided 5 networks each of size 10 and 100 genes[44]. Timeseries data for each network are produced by artificially perturbing a portion of the genes in the network and simulating the response of the network over time.
Specifically, the size 10 DREAM data consist of 10 genes, 21 time points and 5 replicates. The size 100 DREAM data consist of 100 genes, 21 time points, and 10 replicates.
Since our focus is on timeseries data, we did not use the other data sources provided by the DREAM4 challenge. In particular, the DREAM4 challenge provided the results of simulated gene knockout experiments for all genes, and in fact the winning entry in the competition used only the knockout data and ignored the time series data[45]. In practice, however, this is unrealistic, since it is typically feasible to do knockout experiments for only a small proportion of the genes. Time series data do have the potential to provide some information about all potential regulatory relationships, and our goal is to develop methods for doing this, so we have ignored the knockout data here.
Competing methods
To evaluate the performance of our method, we compared it with LASSO[9], as well as with the mutual information methods MRNET, CLR and ARACNE. LASSO is a variable selection method based on a linear regression model, and thus can be used as an alternative to BMA. It is known for being computationally efficient.
We used the implementation of LASSO in the glmnet package in R[12], with the shrinkage penalty parameter, sometimes denoted by λ, chosen by crossvalidation.
Mutual information methods have been used extensively in identifying relationships among genes[21–23]. Since mutual information methods are nondirectional, we used a modified version inspired by[46], including the response as the first column of the matrix given to the mutual information method, and the predictors as the other columns. We then took the first column from the resulting mutual information matrix as the measure of the directed relationship from the predictors to the target gene. MRNET, CLR and ARACNE are different methods for using the mutual information to infer a weighted adjacency matrix between genes. All are implemented in the minet package in R[47].
We have also investigated the possibility of using Dynamic Bayesian Networks (DBN). Methods based on DBNs have been implemented in several R packages. These include the GeneNet package[48], but this is not designed for multiple timeseries and so was not a comparable method on our datasets. The ebdbnet (Empirical Bayes Estimation of Dynamic Bayesian Networks) R package[30], and the G1DBN R package[49], do allow for multiple timeseries, but they are not designed to handle networks of the size of the yeast data. They are more appropriate for networks of sizes up to a few hundred genes. We were able to use the ebdbnet package for the much smaller DREAM4 networks.
YEASTRACT: Assessment of yeast data
The YEASTRACT database[29] is a literaturecurated repository of regulatory relationships between known transcription factors and target genes in yeast, based on more than 1,300 literature references. Among the 3,556 genes in our filtered yeast timeseries data, this database contains documented regulatory relationships for 127 transcription factors (TFs), encompassing a total of 17,173 edges. We therefore evaluate only inferred relationships from these 127 transcription factors, and inferred regulatory relationships from other genes are not used in evaluation.
Availability of supporting data
The yeast time series data are publicly available from ArrayExpresshttp://www.ebi.ac.uk/arrayexpress[50] with accession number EMTAB412. The DREAM4 data are publicly available fromhttp://wiki.c2b2.columbia.edu/dream/index.php?title=D4c2[44].
Abbreviations
 ARACNE:

Algorithm for the reconstruction of accurate cellular networks
 BMA:

Bayesian model averaging
 CLR:

Context likelihood or relatedness
 DBN:

Dynamic Bayesian network
 DREAM:

Dialogue for reverse engineering assessments and methods
 LASSO:

Least absolute shrinkage and selection operator
 MRNET:

Maximum relevance/minimum redundancy
 PRC:

Precision recall curve
 ROC:

Receiver operating characteristic
 TF:

Transcription factor
 YEASTRACT:

Yeast search for transcriptional regulators And consensus tracking.
References
 1.
Murphy K, Mian S: Modelling gene expression data using dynamic Bayesian networks. Technical report, Computer Science Division, University of California, Berkeley, CA, 1999
 2.
Kim S.Y, Imoto S, Miyano S: Inferring gene networks from time series microarray data using dynamic Bayesian networks. Brief Bioinform. 2003, 4 (3): 228235. 10.1093/bib/4.3.228.
 3.
Kim S.Y, Imoto S, Miyano S: Dynamic Bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data. Biosystems. 2004, 75 (1): 5765.
 4.
Zou M, Conzen SD: A new dynamic Bayesian network (DBN) approach for identifying gene regulatory networks from time course microarray data. Bioinformatics. 2005, 21 (1): 7179. 10.1093/bioinformatics/bth463.
 5.
Zhu J, Chen Y, Leonardson AS, Wang K, Lamb JR, Emilsson V, Schadt EE: Characterizing dynamic changes in the human blood transcriptional network. PLOS Comput Biol. 2010, 6 (2): 100067110.1371/journal.pcbi.1000671.
 6.
D’haeseleer P, Wen X, Fuhrman S, Somogyi R: Linear modeling of mRNA expression levels during CNS development and injury. Pacific Symposium on Biocomputing (PSB) conference: January 49, 1999; Hawaii, Volume 4. 1999, 4152.
 7.
Bansal M, Della Gatta G, Di Bernardo D: Inference of gene regulatory networks and compound mode of action from time course gene expression profiles. Bioinformatics. 2006, 22 (7): 815822. 10.1093/bioinformatics/btl003.
 8.
Li Z, Li P, Krishnan A, Liu J: Largescale dynamic gene regulatory network inference combining differential equation models with local dynamic bayesian network analysis. Bioinformatics. 2011, 27 (19): 26862691. 10.1093/bioinformatics/btr454.
 9.
Tibshirani R: Regression shrinkage and selection via the lasso. J R Stat Soc Series B (Methodol). 1996, 58: 267288.
 10.
Tibshirani R, Saunders M, Rosset S, Zhu J, Knight K: Sparsity and smoothness via the fused lasso. J R Stat Soci: Series B (Stat Methodol). 2005, 67 (1): 91108. 10.1111/j.14679868.2005.00490.x.
 11.
Zou H, Hastie T: Regularization and variable selection via the elastic net. J R Stat Soc: Series B (Stat Methodol). 2005, 67 (2): 301320. 10.1111/j.14679868.2005.00503.x.
 12.
Friedman J, Hastie T, Tibshirani R: Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010, 33 (1): 122.
 13.
Raftery AE, Madigan D, Hoeting JA: Bayesian model averaging for linear regression models. J Am Stat Assoc. 1997, 92 (437): 179191. 10.1080/01621459.1997.10473615.
 14.
Hoeting JA, Madigan D, Raftery AE, Volinsky CT: Bayesian model averaging: a tutorial. Stat Sci. 1999, 14: 382401. 10.1214/ss/1009212519.
 15.
Gustafsson M, Hörnquist M, Lundström J, Björkegren J, Tegnér J: Reverse engineering of gene networks with LASSO and nonlinear basis functions. Ann N Y Acad Sci. 2009, 1158 (1): 265275. 10.1111/j.17496632.2008.03764.x.
 16.
Menéndez P, Kourmpetis Y, ter Braak C, van Eeuwijk FA: Gene regulatory networks from multifactorial perturbations using Graphical Lasso: application to the DREAM4 challenge. PLOS ONE. 2010, 5 (12): 1414710.1371/journal.pone.0014147.
 17.
Shojaie A, Michailidis G: Discovering graphical Granger causality using the truncating lasso penalty. Bioinformatics. 2010, 26 (18): 517523. 10.1093/bioinformatics/btq377.
 18.
Yeung KY, Dombek KM, Lo K, Mittler JE, Zhu J, Schadt EE, Bumgarner RE, Raftery AE: Construction of regulatory networks using expression timeseries data of a genotyped population. Proc Nat Acad Sci. 2011, 108 (48): 1943619441. 10.1073/pnas.1116442108.
 19.
Lo K, Raftery A, Dombek K, Zhu J, Schadt E, Bumgarner R, Yeung KY: Integrating external biological knowledge in the construction of regulatory networks from timeseries expression data. BMC Syst Biol. 2012, 6 (1): 10110.1186/175205096101.
 20.
Morrissey ER, Juárez MA, Denby KJ, Burroughs NJ: On reverse engineering of gene interaction networks using time course data with repeated measurements. Bioinformatics. 2010, 26 (18): 23052312. 10.1093/bioinformatics/btq421.
 21.
Basso K, Margolin AA, Stolovitzky G, Klein U, DallaFavera R, Califano A: Reverse engineering of regulatory networks in human b cells. Nat Genet. 2005, 37 (4): 382390. 10.1038/ng1532.
 22.
Margolin AA, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Favera RD, Califano A: ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinformatics. 2006, 7 (Suppl 1): 710.1186/147121057S1S7.
 23.
Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins JJ, Gardner TS: Largescale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLOS Biology. 2007, 5 (1): 810.1371/journal.pbio.0050008.
 24.
Meyer PE, Kontos K, Lafitte F, Bontempi G: Informationtheoretic inference of large transcriptional regulatory networks. EURASIP J Bioinform Syst Biol. 2007, 2007: 79879
 25.
Zoppoli P, Morganella S, Ceccarelli M: TimedelayARACNE: Reverse engineering of gene networks from timecourse data by an information theoretic approach. BMC Bioinformatics. 2010, 11 (1): 15410.1186/1471210511154.
 26.
Lopes F, de Oliveira E, Cesar R: Inference of gene regulatory networks from time series by Tsallis entropy. BMC Syst Biol. 2011, 5 (1): 6110.1186/17520509561.
 27.
Zellner A: On assessing prior distributions and Bayesian regression analysis with gprior distributions. Bayesian Inference Decis Tech: Essays Honor of Bruno De Finetti. 1986, 6: 233243.
 28.
Madigan D, Raftery AE: Model selection and accounting for model uncertainty in graphical models using Occam’s window. J Am Stat Assoc. 1994, 89 (428): 15351546. 10.1080/01621459.1994.10476894.
 29.
Teixeira MC, Monteiro P, Jain P, Tenreiro S, Fernandes AR, Mira NP, Alenquer M, Freitas AT, Oliveira AL, SáCorreia I: The YEASTRACT database: a tool for the analysis of transcription regulatory associations in saccharomyces cerevisiae. Nucleic Acids Res. 2006, 34 (suppl 1): 446451.
 30.
Rau A, Jaffrézic F, Foulley J, Doerge RW: An empirical Bayesian method for estimating biological networks from temporal microarray data. Stat Appl Genet Mol Biol. 2010, 9 (1): 15446115.
 31.
Bottolo L, Richardson S: Evolutionary stochastic search for Bayesian model exploration. Bayesian Anal. 2010, 5 (3): 583618. 10.1214/10BA523.
 32.
Yeung KY, Bumgarner RE, Raftery AE: Bayesian model averaging: development of an improved multiclass, gene selection and classification tool for microarray data. Bioinformatics. 2005, 21 (10): 23942402. 10.1093/bioinformatics/bti319.
 33.
Guelzim N, Bottani S, Bourgine P, Képès F: Topological and causal structure of the yeast transcriptional regulatory network. Nat Genet. 2002, 31 (1): 6063. 10.1038/ng873.
 34.
Zhu J, Zhang B, Smith EN, Drees B, Brem RB, Kruglyak L, Bumgarner RE, Schadt EE: Integrating largescale functional genomic data to dissect the complexity of yeast regulatory networks. Nat Genet. 2008, 40 (7): 854861. 10.1038/ng.167.
 35.
Yip KY, Alexander RP, Yan K, Gerstein M: Improved reconstruction of in silico gene regulatory networks by integrating knockout and perturbation data. PLOS ONE. 2010, 5 (1): 812110.1371/journal.pone.0008121.
 36.
Clyde M, George E. I: Model uncertainty. Stat Sci. 2004, 19: 8194. 10.1214/088342304000000035.
 37.
Raftery AE: Bayesian model selection in social research. Sociol Methodol. 1995, 25: 111164.
 38.
Raftery AE: Bayes factors and BIC. Sociol Methods Res. 1999, 27 (3): 411417. 10.1177/0049124199027003005.
 39.
Dempster AP, Laird NM, Rubin DB: Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc Series B (Methodological). 1977, 39 (1): 138.
 40.
McLachlan G, Krishnan T: The EM Algorithm and Extensions. Volume 382. 2007, Hoboken, New Jersey: John Wiley & Sons
 41.
Marbach D, Prill RJ, Schaffter T, Mattiussi C, Floreano D, Stolovitzky G: Revealing strengths and weaknesses of methods for gene network inference. Proc Nat Acad Sci. 2010, 107 (14): 62866291. 10.1073/pnas.0913357107.
 42.
Marbach D, Schaffter T, Mattiussi C, Floreano D: Generating realistic in silico gene networks for performance assessment of reverse engineering methods. J Comput Biol. 2009, 16 (2): 229239. 10.1089/cmb.2008.09TT.
 43.
Prill RJ, Marbach D, SaezRodriguez J, Sorger PK, Alexopoulos LG, Xue X, Clarke ND, AltanBonnet G, Stolovitzky G: Towards a rigorous assessment of systems biology models: the DREAM3 challenges. PLOS ONE. 2010, 5 (2): 920210.1371/journal.pone.0009202.
 44.
DREAM4 In Silico Network Challenge. [http://wiki.c2b2.columbia.edu/dream/index.php?title=D4c2],
 45.
Pinna A, Soranzo N, de la Fuente A: From knockouts to networks: establishing direct causeeffect relationships through graph analysis. PLOS ONE. 2010, 5 (10): 1291210.1371/journal.pone.0012912.
 46.
Shimamura T, Imoto S, Yamaguchi R, Fujita A, Nagasaki M, Miyano S: Recursive regularization for inferring gene networks from timecourse gene expression profiles. BMC Syst Biol. 2009, 3 (1): 4110.1186/17520509341.
 47.
Meyer PE, Lafitte F, Bontempi G: minet: AR/Bioconductor package for inferring large transcriptional networks using mutual information. BMC Bioinformatics. 2008, 9 (1): 46110.1186/147121059461.
 48.
Schäfer J, OpgenRhein R, Strimmer K: Reverse engineering genetic networks using the Genenet package. J Am Stat Assoc. 2001, 96: 11511160. 10.1198/016214501753382129.
 49.
Smith SM, Fulton DC, Chia T, Thorneycroft D, Chapple A, Dunstan H, Hylton C, Zeeman SC, Smith AM: Diurnal changes in the transcriptome encoding enzymes of starch metabolism provide evidence for both transcriptional and posttranscriptional regulation of starch metabolism in Arabidopsis leaves. Plant Physiol. 2004, 136 (1): 26872699. 10.1104/pp.104.044347.
 50.
ArrayExpress. [http://www.ebi.ac.uk/arrayexpress],
Acknowledgements
We would like to thank Drs. Roger Bumgarner and Chris Fraley for helpful discussion.
Funding
AER, KYY and WCY are supported by NIH grant 5R01GM084163. KYY is also supported by 3R01GM08416305S1. AER is also supported by NIH grants R01 HD054511 and R01 HD070936, and by Science Foundation Ireland ETS Walton visitor award 11/W.1/I207.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
WCY participated in method development, implemented the ScanBMA method, carried out empirical studies comparing ScanBMA to other competing methods, and drafted the manuscript. AER designed the ScanBMA algorithm, designed the empirical studies and edited the manuscript. KYY identified the datasets, participated in the design of the empirical studies, and assisted in manuscript preparation. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Bayesian inference
 Bayesian model averaging
 Gene regulatory networks