key: cord-0015322-4ppst4dg
authors: Lu, Jonathan; Dumitrascu, Bianca; McDowell, Ian C.; Jo, Brian; Barrera, Alejandro; Hong, Linda K.; Leichter, Sarah M.; Reddy, Timothy E.; Engelhardt, Barbara E.
title: Causal network inference from gene transcriptional time-series response to glucocorticoids
date: 2021-01-29
journal: PLoS Comput Biol
DOI: 10.1371/journal.pcbi.1008223
sha: 3936ba471cd2cb26e9634d1a115a376e21225ddd
doc_id: 15322
cord_uid: 4ppst4dg

Gene regulatory network inference is essential to uncover complex relationships among gene pathways and inform downstream experiments, ultimately enabling regulatory network re-engineering. Network inference from transcriptional time-series data requires accurate, interpretable, and efficient determination of causal relationships among thousands of genes. Here, we develop Bootstrap Elastic net regression from Time Series (BETS), a statistical framework based on Granger causality for the recovery of a directed gene network from transcriptional time-series data. BETS uses elastic net regression and stability selection from bootstrapped samples to infer causal relationships among genes. BETS is highly parallelized, enabling efficient analysis of large transcriptional data sets. We show competitive accuracy on a community benchmark, the DREAM4 100-gene network inference challenge, where BETS is one of the fastest among methods of similar performance and additionally infers whether causal effects are activating or inhibitory. We apply BETS to transcriptional time-series data of differentially-expressed genes from A549 cells exposed to glucocorticoids over a period of 12 hours. We identify a network of 2768 genes and 31,945 directed edges (FDR ≤ 0.2). We validate inferred causal network edges using two external data sources: Overexpression experiments on the same glucocorticoid system, and genetic variants associated with inferred edges in primary lung tissue in the Genotype-Tissue Expression (GTEx) v6 project. BETS is available as an open source software package at https://github.com/lujonathanh/BETS.

The recent availability of gene expression measurements over time has enabled the search for interpretable statistical models of gene regulatory dynamics [1] . These time-series data present a unique opportunity to use the coordinated transcriptional response to environmental exposure to infer causal relationships between genes. However, there are several challenges to overcome in the analysis of time-series transcriptomic data. These data are generally highdimensional: the number of quantified gene transcripts-approximately 20,000 in human samples-often dramatically exceeds the number of available time points and samples. Many classical statistical assumptions fail to hold in this high-dimensional regime [2, 3] . Moreover, the large number of gene transcripts poses a computational burden, as the number of possible edges in a gene network grows quadratically. Also, a transcriptional time series often has a small number of time points, and those time points are often not uniformly spaced; furthermore, because transcriptional time-series data often quantify transcription post exposure, the time series is not stationary, and genes respond to the exposure and return to baseline at different rates [4, 5] .

In this work, we develop an approach that uses gene transcription time series following glucocorticoid (GC) exposure to build a directed gene network [6] . GCs play an essential role in regulating stress response, and are widely used as anti-inflammatory and immunosuppresive medications [6, 7] . Dexamethasone and other GCs have recently been recommended by the U.K. Government, U.S. National Institutes of Health, World Health Organization, and Infectious Disease Society of America for treatment of hospitalized COVID-19 patients with severe disease who are on mechanical ventilation or extracorporeal membrane oxygenation [8] [9] [10] [11] . Despite clinical benefits, prolonged exposure to GCs has been linked to increased risk of type 2 diabetes mellitus (T2DM) [12] and obesity [13] . Understanding the immune, metabolic, and transcriptional effects of GCs may enable the development of improved anti-inflammatory treatments without metabolic side effects. A recent study assayed A549 lung cells over 12 hours to characterize the effect of GCs on cell state [6] . Here, we develop a method to accurately, interpretably, and efficiently infer a directed gene network using the study's transcription time-series data. We focus our network analysis on immune-related genes, metabolism-related genes, and transcription factors (TFs) to study the inferred coordinated response of these systems to GCs.

Our method, Bootstrap Elastic net inference from Time Series (BETS), uses vector autoregression with elastic net regularization to infer directed edges between genes. Stability selection, which assesses the robustness of an edge to perturbations in the data, leads to improvements over baseline vector autoregression methods in this high-dimensional context [3] . Furthermore, BETS is biologically interpretable because estimated coefficients provide the direction (sign) and effect size of the causal relationship between a pair of genes. Finally, BETS's parallelization enables efficient inference of networks with millions of possible edges in a computationally tractable way.

We use the causal network inferred by BETS on the GC time-series data to study the relationships between TFs, immune genes, and metabolic genes. We validate our network using two approaches: Ten measurements of the same GC system with an overexpressed TF, and an expression quantitative trait loci (eQTL) study in human primary lung tissue [14] . Although our framework is motivated by transcriptional response to GC exposure, our approaches are general, and BETS is able to infer directed networks from arbitrary high-dimensional timeseries data.

Several methods have been developed to estimate directed gene networks from transcription time-series data (S1 Fig) [15] [16] [17] [18] [19] [20] [21] [22] [23] . These methods estimate directed networks in which the directed edges between nodes-representing genes-indicate a cause-effect relationship between genes. In other words, perturbing expression of the causal gene would lead to changes in expression of the effect gene [24] . We briefly overview these methods; for detailed discussion, see the Supplemental Information. Here, we take g 0 to be the causal gene and g to be the effect gene, and we quantify support for a causal edge g 0 ! g in the time-series data.

Mutual information (MI) methods assess the MI between the expression of g 0 at the previous time point and the expression of g at the current time point (S1(A) Fig) [25] [26] [27] [28] [29] [30] . A causal edge g 0 ! g is included in the network if the MI of the two genes across time exceeds a threshold.

Granger causality methods determine if including the expression of g 0 at the previous time point improves our ability to predict the expression of g at the current time point beyond using the expression of g at the previous time point [31] . A common way to implement Granger causality is through a vector autoregression (VAR) model, which assumes a linear relationship between all genes' expression at the previous time point and the expression of g at the current time point. A causal edge g 0 ! g is included in the network when g 0 has a statistically significant coefficient in the VAR.

Ordinary differential equations (ODEs) fit the derivative of the expression of g as a function of all genes' expression at a single time point (S1(C) Fig) [15, 32, 33] . ODE methods typically assume linearity, as small sample sizes make it challenging to infer the parameters of nonlinear functions. A causal edge g 0 ! g is included when g 0 has a statistically significant coefficient in the ODE.

Decision trees (DTs) are a type of nonparametric function based on partitioning the data [34, 35] . DT methods fall either under VAR or ODE; either the DTs fit the expression of g at the current time as a function of all genes' expression at the previous time point (VAR), or they fit the derivative of the expression of g as a function of all genes' expression at a single time point (ODE) (S1(D) Fig) [36, 37] . A causal edge g 0 ! g is included in the network when an importance score for g 0 exceeds some threshold, where importance scores are typically the reduction in variance of g when g 0 is included as a predictor.

Dynamic Bayesian networks (DBNs) search the space of possible directed acyclic graphs between previous and current expression levels to identify the network structure with the highest posterior probability of each edge given the data (S1(E) Fig) [38] [39] [40] [41] [42] . DBNs typically assume a linear relationship between previous and current expression. A causal edge g 0 ! g is included in the network when its marginal posterior probability of existence exceeds some threshold.

A Gaussian process (GP) is a distribution over continuous, nonlinear functions. GPs are often used in the context of nonlinear DBNs, where GP regression is used to model a nonlinear relationship between previous expression and current expression (S1(F) Fig) [43, 44] . A causal edge g 0 ! g is included in the network when its posterior probability of existence exceeds some threshold.

While these approaches produce directed networks that have the flavor of Bayesian networks, except for DBNs, none of them produce graphs that are constrained to be acyclic, so they do not have the same statistical semantics as Bayesian networks.

First, we briefly describe the approach that BETS uses to infer a directed gene network. Next, we compare results from BETS to those from twenty two other methods on the 100-gene timeseries data from the DREAM4 Network Inference Challenge [45] . Then, we describe the network estimated from the GC transcription time-series data. Finally, we validate the inferred network using two different frameworks: Overexpression experiments on the same system, and genetic variants associated with inferred edges in human primary lung tissue in the Genotype-Tissue Expression (GTEx) v6 project [14] .

Directed networks represent causal relationships among diverse interacting variables in complex systems. We developed a robust, scalable approach based on ideas from Granger causality to construct these directed networks from short, high-dimensional time-series observations of gene expression levels.

Let G be the set of all p = |G| genes in the data set and g 2 G be a gene. Let ¬g be G with g removed. Let t be a single time point, ranging from {1, 2, . . ., T}. Let X g t be the expression of gene g at time t. Let L be the time lag, or the number of previous time point observations; so L = 2 means that we use two previous time points, t − 1 and t − 2, to predict expression at time t. These types of autoregressive models work best with similarly-spaced time points, as the data sets in this paper approximate, and assume stationarity, or the same causal effects across each time gap. Definition 2.1 (Granger causality) . For lag L, a gene g 0 is said to Granger-cause another gene g if using X g 0 tÀ 1 ; . . . ; X g 0 tÀ L , the expression values of g 0 at times t − 1 to t − L, improves prediction of X g t , the expression value of g at time t, beyond predicting X g t using X g tÀ 1 ; . . . ; X g tÀ L alone.

To test for Granger causality from g 0 to g, we first preprocessed the gene expression timeseries data (Methods). For every potential effect gene g, we fit all other genes g 0 2 ¬g simultaneously (Eq 1), echoing ideas from the graphical lasso for undirected network inference [46] . Intuitively, this adapts the idea of Granger causality to conditional Granger causality, where we consider how gene g 0 Granger-causes g conditioning on the effects of all other genes. This approach uses the regression:

where � t � N ð0; 1Þ. For BETS, we set L = 2. To test for an edge, if there is statistical support for b g 0 ;g ' 6 ¼ 0, then we say g 0 conditionally Granger-causes g at lag L. We build the directed network by including a directed edge to g from every gene g 0 that has been inferred to conditionally Granger-cause g.

Robustly building this network is difficult due to the high dimensionality of the problem: The number of genes that could Granger-cause a given g far exceeds the available time points and technical replicates. To address this challenge, BETS regularizes the VAR model parameters using an elastic net penalty (Methods, Fig 1A) . Elastic net regression encourages sparsity and performs automatic variable selection on the genes being tested for causal influence [47] . The elastic net penalty, unlike the lasso penalty [48] , is able to select groups of correlated variables and allows the number of selected variables to be greater than the number of samples. This is important for gene expression assays where gene expression levels are often well correlated, and there are far more genes than samples [2] .

In BETS, we fit the same VAR model to a data set in which causal genes have their expression permuted over time to generate a null distribution of edge coefficients. The coefficients are thresholded to produce a causal network with each edge at edge false discovery rate (FDR) � 0.05 ( Fig 1A) . We then applied this network inference procedure to multiple (here, 1000) bootstrapped samples of the original data set ( Fig 1B) . Each edge has a selection frequency, or the frequency that the edge appears in networks inferred from the bootstrapped samples. Inspired by stability selection, this approach assesses if network edges are robust to perturbations of the data [3] . Finally, we ran this overall procedure on a permuted version of the original data set to obtain a null distribution of selection frequencies (Fig 1C) . The selection frequency threshold for including each edge is chosen to control the stability FDR �0.2. As a baseline, we compare BETS against Enet, which runs elastic net regression without stability selection to produce a causal network with each edge at edge FDR � 0.05 (Fig 1A) .

We evaluated BETS against other directed network inference methods. We used the DREAM4 Network Inference Challenge [45] , a community benchmark for directed network inference using gene time-series data. This benchmark consisted of five data sets, each with ten timeseries measurements for 100 genes across 21 time points [45] . Evaluation was previously done by looking at the average of the area under the precision recall curve (AUPR) or the area under the receiver operating characteristic (AUROC) over the five data sets [37, 45] . Any method that provides a ranking of possible network edges could be evaluated in this framework.

We tested BETS and Enet against 22 other methods on the DREAM challenge [36, 37, 40, [49] [50] [51] . We ran SWING-RF, SWING-Lasso, CSId, Jump3, CLR, MRNET, and ARACNE inhouse and found our results consistent with those reported in the literature. All 22 methods reported AUPR, but only 17 reported AUROC.

BETS ranked 7th out of 24 in AUPR with an average AUPR of 0.128 (Fig 2A and S1 Table) and 4th out of 19 in AUROC with an average AUROC of 0.688 (Fig 2B and S2 Table) . BETS was the top performer of all VAR methods, and Enet was second best. All 24 methods outperformed random selection of edges, which achieved an average AUPR of 0.002 and average AUROC of 0.50 [49] . We also found that BETS and Enet had similar performance to the DBN methods in AUPR, and outperformed most of them in AUROC. Ranked by the top AUPR of each class of methods, the best performing class was DT, followed by GP, MI, VAR, DBN, and ODE [36, 40, 49] . The VAR method used in BETS produces edge signs (indicating excitatory or inhibitory causal effects) and effect sizes. While other methods based on GPs (e.g., CSId), MI (e.g., tl-CLR) or DTs (e.g., SWING-RF) had marginally better overall network inference, they do not provide insight into the causal relationships because they only output a positive measure of a causal interaction [32, 44, 51] .

Next, we compared the speed of BETS and three other top-performing methods: SWING-RF, CSId and Jump3 (S3 Table) . SWING-RF was the fastest at 0.11 hours, while BETS took 4.8 hours, CSId took 9.8 hours and Jump3 took 45 hours. Thus, while BETS had a lower AUPR compared to CSId and Jump3, it was substantially faster. BETS had both a lower AUPR and longer runtime than SWING-RF.

BETS improved upon Enet using stability selection. To quantify this improvement, we compared three other models: Elastic net with lag L = 1, ridge regression with lag L = 2, and lasso with lag L = 2 (S4 Table) . In each case, the stability selection version outperformed the original version in average AUPR and AUROC. The improvement in average AUPR ranged between 0.016 and 0.03 (+20% to +31%), while the improvement in average AUROC ranged between 0.012 and 0.04 (+1.8% to +6.1%). Thus, our stability selection procedure leads to improved performance for multiple versions of VAR. We also found that stability selection performance is robust to the number of bootstrap samples (S5 Table) . Decreasing the number of bootstrap samples from 1000 to 100 led to minor decreases of −0.004 in AUPR and −0.008 in AUROC, within the standard deviation across the networks. It also resulted in a 10-fold decrease in memory usage and 3-fold decrease in run time, due to a constant-time hyperparameter search. If users face computational constraints, we recommend that they use 100 bootstrap samples for nearly equivalent performance.

Finally, we found that BETS' performance on DREAM is robust to the choice of lag. We ran BETS with lag L = 1 instead of lag L = 2 (S4 Table) . BETS with lag L = 1 achieves an average AUPR of 0.14 (an increase of 0.012 from lag L = 2) and an average AUROC of 0.686 (a decrease of 0.002 from lag L = 2, within the standard deviation across networks). BETS with lag L = 1 still ranks 7th out of 24 in AUPR and 4th out of 19 in AUROC (tied with GP4GRN). Thus, BETS achieves consistently good performance on DREAM for both lags L = 1 and L = 2.

To infer the causal relationships in the GC response network, we analyzed RNA-seq data collected from human adenocarcinoma and lung model cell line A549, which consists of two data sets. In an original exposure data set, cells were exposed to the synthetic GC dexamethasone (dex) for 0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 10, and 12 hours [6] . In an unperturbed data set, the cells were first exposed to dex for 12 hours, after which the media was replaced and dex removed, and then measurements were taken at the same intervals 0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 10, and 12 hours. BETS was fit jointly over the two data sets. In total there were 7 technical replicates (4 from original exposure and 3 from unperturbed). A single VAR was fit on 70 samples: Each of the 7 replicates had 10 samples, because using a lag L = 2 VAR model turns 12 time points into 10 samples.

We applied BETS to the GC-mediated expression responses to infer a causal network ( Fig  3A) . Edges with selection frequency (i.e., frequency of appearance among bootstrap networks) at least 0.097 were declared significant (FDR � 0.2; Fig 3B) . The network contained 2, 768 nodes representing distinct genes and 31, 945 directed edges (0.4% of possible edges). Of these, 466 genes were causes (i.e., had an outward directed edge) and all 2, 768 genes were effects (i.e., had an incoming directed edge). In Granger causality, and dynamical systems more generally, a causal gene g 0 is allowed to have incoming directed edges because g 0 may be affected by the past value of another gene g@, and g 0 may have a causal effect on the later value of gene g. The out-degree distribution was heavy-tailed and skewed right (Fig 3C) , while the in-degree distribution was lighter-tailed and more symmetric (Fig 3D) . The network's edge in-degree had a heavier left tail and lighter right tail than a normal distribution ( Fig 3E) . This suggests that causal genes are relatively rare (only 1/6th of network genes are causes) and a fifth of those only affect a single gene, whereas genes that are effects tend to have multiple causes. The network was inferred efficiently due to parallelization across genes, taking six days in real time and 292 days in CPU time to fit 5.5 million elastic net models.

To study the network with respect to the glucocorticoid system, we annotated specific genes as transcription factors (TFs), immune-related, or metabolism-related [52] [53] [54] [55] . First, we inspected enrichment of each category among the causal genes ( Fig 3F) . At FDR � 0.05, we found enrichment for TFs among causes; there were 226 causal TFs, representing 8.2% of the 2, 768 input genes. 62 of these TFs were causal, representing 13% of all causal genes (odds ratio (OR) = 2.0, Fisher's exact test (FET) adjusted p � 2.9 × 10 −5 ). Similarly, we found an enrichment among immune-related genes as causes: of 109 immune genes, representing 3.9% of the input genes, 39 of these were causes, representing 8.4% of all causal Table. https://doi.org/10.1371/journal.pcbi.1008223.g003 genes (OR = 2.9, FET adjusted p � 2.5 × 10 −6 ). In contrast, there was no enrichment among metabolism-related genes: there were 120 metabolic genes, representing 4.3% of input genes; 19 of these metabolism genes were causes, representing 4.1% of all causes (OR = 0.93, FET adjusted p � 0.66). This suggests that our network is enriched for causal TFs and immune genes.

To study the interactions among gene classes inferred by our network, we quantified enrichment for edges between each of the four gene classes-immune, metabolic, TF, and other gene types (Fig 3G; S6 Table) . We found enrichment of 11 of the 16 possible edge types (FDR � 0.05). The network was enriched for edges from i) causal TFs to immune genes and other genes; ii) causal immune genes to TFs, immune genes, metabolic genes, and other genes; and iii) causal metabolic genes to TFs, immune genes, metabolic genes, and other genes. This suggests that our network is enriched for a broad range of edge types.

The inferred relationships in BETS are conditionally Granger-causal, meaning that we find that gene g 0 Granger-causes g conditioning on the effects of other genes. Thus, an edge g 0 ! g should be assessed based on g's residuals after controlling for the effects of other genes and itself, instead of g's raw values. Consider the edge KRT6A ! NKAIN4 (Fig 4) : NKAIN4's raw expression values suggest a negative relationship between KRT6A and NKAIN4 (Fig 4A and  4B ). However, after controlling for the effects of all covariates besides KRT6A, a positive relationship between KRT6A and NKAIN4 appears (Fig 4C and 4D) . Thus, conditionally Grangercausal relationships g 0 ! g should assess g only after controlling all other effects on g.

Our network identified known biological interactions between genes with immune, metabolic, and TF roles; we highlighted 16 gene pairs with experimentally validated interactions ( Fig 5, S2 Fig, and S7 Table) . Known interactions were found using the BIOGRID PPI database [56] . SOCS1 and SOCS3 bind IRS2 and promote its degradation, leading to reduced insulin signaling [57, 58] ; furthermore, SOCS1 represses IL-4-induced IRS2 singling [59] . NR4A1 heterodimerizes with RXRA to activate it in order to promote gene expression under vitamin A signaling [60] ; NR4A1 also inhibits p300-induced RXRA acetylation [61] . Eleven of the 16 edges had the correct interaction direction; the five that were reversed are TNFAIP3 ! IRAK2, SOCS3 ! HIVEP1, ATF3 ! MDM2, E2F1 ! CDH1, and FOS ! EGFR. These results suggest that BETS infers biologically meaningful relationships, but transcriptional data, absent other assays on protein abundance and cellular dynamics, are often underpowered to resolve the direction of the edge.

We asked whether our inferred network edges validated on overexpression versions of the same experimental system, in which each of ten TFs was separately overexpressed over the same 12 hours of observations. Specifically, we assessed the concordance between inferred network edges g 0 ! g and their coefficient in the overexpression data set under a VAR model (Methods).

We first evaluated how well network edges replicated on individual overexpression data sets. We performed linear regression of a one-hot encoding of the original network's edge sign (i.e., positive versus no edge or negative; negative versus positive or no edge) as the predictor against the VAR model edge coefficients estimated from each of the overexpression time series as the response (Fig 6A and 6B , Methods). Of the ten data sets, 9 showed enriched positive effect sizes among positive edges at FDR �0.2 (CEBPB, CEBPD, FOSL2, FOXO1, FOXO3, KLF6, KLF9, KLF15, OCT4; Fig 6A) . Three data sets showed enriched negative effect sizes among negative edges (OCT4, TFCP2L1, CEBPD) and four showed enriched positive effect sizes among negative edges (CEBPB, FOSL2, KLF9, KLF15; Fig 6B) .

Taken together, the positive edges inferred by BETS validate on the overexpression data, but the negative edges do not, indicating repressive effects may have inconsistent signs or feedback loops.

Next, we checked whether the 123 inferred causal edges from the TF TFCP2L1 validated in the TFCP2L1 overexpression data set (there were only about 10 causal edges from each of the other 9 TFs). We regressed the original network's edge sign (+ 1 for positive edges, 0 for no edge, and −1 for negative edge) as the predictor against the overexpression VAR model edge coefficients as the response (Fig 6C) . We found a positive relationship between the edge sign and overexpression coefficient (slope 0.17, two-sided t-test p � 5 × 10 −5 ). This shows that causal edges from TFCP2L1 are enriched for matched effect directions in the TFCP2L1 overexpression data.

This validation may be limited by a misspecification of the linear regression model. As a supplementary analysis, we fit nonstationary GPs to gene trajectories using nsgp; genes that 

Causal gene network inference from time series were differentially expressed in the TF overexpression data compared to the original exposure data were listed as potential targets of that TF (Supplemental Information). We found that BETS weakly predicts potential targets of 5 of the 10 over-expression TFs. This approach did not show substantial improvement above random prediction of potential targets on these data (FDR � 0.1).

We next validated our network edges on an expression quantitative trait-loci (eQTL) study. A single nucleotide polymorphism (SNP) S is an eQTL for a gene g 0 if it is associated with g 0 's expression level within a population. Given a true causal edge g 0 ! g, if a SNP S is a local (cis-) eQTL for g 0 , S might also be a distal (trans-) eQTL for g [62] . We used gene expression levels in primary lung tissue (n = 278) from the Genotype Tissue Expression (GTEx) project v6p [14] . We observed an enrichment of low trans-eQTL association p-values from the directed network compared to shuffling the variant labels (Fig 7A and 7B ). This suggests our network captures more valid causal effects than expected by chance.

We next inspected specific associations and their corresponding edges. We found 340 trans-eQTL pairs in lung samples corresponding to 130 network edges (q-value FDR � 0.2). There are more trans-eQTLs than edges because there are multiple cis-eQTLs for some causal genes g 0 . The 340 trans-eQTLs greatly improved upon the 2 identified in primary lung tissue in the GTEx v6 trans-eQTL study [63] , demonstrating the utility of transcriptional time series for prioritizing promising associations. The top trans-associations were rs2302178-CLDN1 (q-value FDR � 0.095; extended from the cis-association rs2302178-HS3ST6), rs590429-ADAMTS 

(q-value FDR � 0.11; extended from the cis-association rs2302178-OLR1), and rs2072783-CLIP2 (q-value FDR � 0.11; extended from the cis-association rs2302178-GMPR; Fig 7C) .

We searched for validated associations between immune-related genes, metabolic-related genes, and TFs. One association was OLR1 ! ITGAV, where the known association between SNP rs4329754 and OLR1 extends to an association between the same genetic variant and effect gene ITGAV (q-value FDR �0.13) [14] . OLR1 plays key roles in immunity and metabolism [64, 65] . It is associated with metabolic syndrome [66] and atherosclerosis [66] , and modulates inflammatory and humoral immune responses [67, 68] . Meanwhile, ITGAV plays a key role in the motility of CD4 + T cells during inflammation [69] .

Another association was between the TF SNAI2 and gene PTPN6, where we find that the known association between genetic variant rs56800165 and SNAI2 extends to an association between the same SNP rs56800165 and effect gene PTPN6 (q-value FDR �0.17) [14] . SNAI2 is a direct target of the glucocorticoid receptor GR that regulates cell migration in breast cancer [70] , while PTPN6 is involved in glucose homeostasis via negative regulation of insulin signalling [71] . PTPN6 is also associated with inflammatory phenotypes in multiple diseases [72, 73] . Finally, both SNAI2 and PTPN6 are involved in the cell-cell adherens junctions pathway, as SNAI2 represses transcription of cadherin, while PTPN6 positively regulates the cadherin-catenin complex [74] . Thus, for several eQTL-validated edges for gene pairs, we find that the genes are involved in related biological processes, but further experimentation is required to confirm direct interactions.

As A549 cells are models for lung tissue [75] , we quantified enrichment of validated edges in lung compared to enrichment in four non-lung tissues with similar sample sizes: subcutaneous adipose (n = 298), transformed fibroblasts (n = 272), tibial artery (n = 285), and thyroid (n = 278). We validated 341 unique network edges across the five tissues (FDR � 0.2). 130 edges validated for lung, 4 for subcutaneous adipose, 125 for transformed fibroblasts, 3 for tibial artery, and 82 for thyroid tissues. More network edges validated in primary lung than in non-lung tissues, suggesting that A549 cells most closely match lung samples among GTEx tissues; this is consistent with their tissue of tumor origin.

We described an approach, BETS, to build directed networks using short time-series observations of high-dimensional transcription data. BETS combined ideas from elastic net regression, graphical lasso, stability selection, and VAR models to infer Granger-causal relationships in high-dimensional transcription time-series data. Our method achieved competitive performance on the DREAM4 100-Gene Network Inference Challenge, ranking 7th out of 24 methods in AUPR and 4th out of 19 methods in AUROC; it was also faster than several methods with similar or better performance and infers effect size and sign, unlike the other top performing methods. Stability selection resulted in consistent improvement to VAR models across different hyperparameter settings.

Next, we applied BETS to time-series RNA-seq data from human A549 cells exposed to glucocorticoids and identified a directed network of 31, 945 edges (FDR �0.2), capturing the causal relationships among genes after exposure to GCs. Despite the intervention of GCs in this cell line, we expect that many of the causal relationships estimated using these data will exist across cellular environments in lung cells. In our estimated causal network, we found enrichment of immune genes and TFs among causal genes. We also found enrichment of 11 specific types of causal edges from TFs, immune genes, and metabolic genes. We validated our network first in ten overexpression data sets. Edges that were positive in the original network had an enrichment of positive VAR effect sizes in the overexpression data. However, edges in the original network did not predict differential expression of genes in the overexpression data, as called by nsgp. Validating network edges by searching for trans-eQTLs in GTEx primary lung tissue samples, we found an enrichment of associations with genetic variants across network edges. Finally, we discovered 340 trans-eQTLs, a dramatic improvement from the GTEx v6 trans-eQTL study [63] .

While BETS has demonstrated effective inference of causal relationships, there are interesting future directions to explore. All methods that infer networks from transcriptional time series face several difficulties: i) Transcript levels are sometimes an imperfect proxy for protein levels, especially when transcript dynamics are changing [15, 76] ; ii) the scarcity of time point samples causes statistical challenges for inferring millions of possible causal interactions between genes, let alone non-additive interactions among causes [3, 77, 78] ; iii) transcription data do not capture the complete regulatory context including chromatin structure and epigenetic regulations [15] ; iv) transcription relationships are often nonstationary: the relationship may change over time due to responses from the environment [4, 5] ; and v) inferred networks are often sensitive to the choice of preprocessing and parameter choices [79] . Single-cell data also implicitly include transcription time-series information when pseudotime is inferred, making ideas from Granger causality exceptionally relevant. However, recent preprints suggest some limits to pseudotime's potential for network inference [80, 81] . Finally, experimental followup is key to establishing causality; BETS only generates promising, interpretable hypotheses. Indeed, by discovering hundreds more trans-eQTLs than the GTEx study (a 170-fold increase) [63] , BETS demonstrates its potential to prioritize biologically meaningful associations.

Bootstrap Elastic net regression from Time Series (BETS) is a vector-autoregressive approach to causal inference from gene expression time-series data. It is based on the principle of Granger causality [31] : a gene g 0 Granger-causes another gene g if previous information from gene g 0 improves our current predictions of gene g, beyond using previous information from g and from other genes.

BETS first preprocesses the data. BETS fits an elastic net vector autoregression model to handle the high dimensionality of the time series, inferring a network (Fig 1A) . It infers one network for each of 1000 bootstrapped samples of the original data set and computes each edge's selection frequency, or its frequency of appearance among the bootstrapped networks ( Fig 1B) [3] . Finally, BETS includes an edge in the network using the selection frequencies ( Fig  1B) . Our baseline comparison, Enet, only preprocesses the data and fits an elastic net vector autoregression model from the original data (Fig 1A; Section 2 

Preprocessing temporal time-series data. For a gene temporal profile (i.e., one gene's expression values across time for a single replicate), we used zero-mean unstandardized normalization, which centers each gene temporal profile to have mean zero across time. Because gene temporal profile ranges from staying almost constant to having drastic fluctuations, BETS uses this approach because a unit-variance normalization would over-represent the weak causal effects of genes with lower variability.

Vector autoregression model. Let G be the set of all genes in the data, let p = |G| be the number of genes, and let g be a gene. Let ¬g be G with g removed. Let there be T time points total, and let t 2 {1, 2, . . ., T} be a single time point. Let there be R replicates of the gene expression time series. Let X g t;r be the expression of gene g at time t for replicate r. Let X g t ¼ ½X g t;1 ; X g t;2 ; . . . ; X g t;R � T be the R × 1 vector of gene expression levels of gene g across R replicates at time t. The rest of the paper does not mention replicates for simplicity, but here we discuss replicates for completeness. Let g 0 be the gene we are testing to be causal for gene g and let ℓ refer to the time lag of the causal edge g 0 ! g. Let L be the maximum lag. In BETS, L = 2 is the default.

We model each gene g as

where � t � N ð0; 1Þ. In other words, the expression of each gene g is modeled as a linear function of its and other genes' L previous expression values, under independent Gaussian noise. a g ' represents the (scalar) effect size of gene g's ℓth previous value, X g tÀ ' , on its current value, X g t . b g 0 ;g ' represents the (scalar) effect size of the ℓth previous value of gene g 0 6 ¼ g, X g 0 tÀ ' , on gene g's current value, X g t . Eq 2 requires that t > ℓ for the ℓth previous value, X g tÀ ' , to exist. To demonstrate how our model is fit in practice, we reformulate Eq 2 using matrix notation. Each row represents one time point for one replicate. There are T − L time points with t > L and R replicates, so there are R(T − L) samples, or rows, in total. Let N = R(T − L).

Define X g t , an N × 1 vector, as: 

We can similarly write X g tÀ ' , which is X g t with each entry replaced by its ℓth previous value. Define X g tÀ ' , a N × L matrix consisting of the first L previous vectors X g tÀ ' , i.e., for ℓ ranging in {1, . . ., L}.

Let a g ' be a L × 1 vector of the L lagged coefficients. 

Next, let us formulate Eq 2 involving the genes g 0 in matrix notation. Let X :g tÀ ' be a N × L(|G| − 1) predictor matrix of the vectors X g 0 tÀ ' , for g 0 tÀ ' .

Let b :;g ' be a L(|G| − 1) × 1 vector of the causal coefficients b g 0 ;g ' where g 0 6 ¼ g.

. . 

We then fit the model:

where � t is a N × 1 vector with each element � t;n � N ð0; 1Þ. In its most compact form, we can write

Note that X G tÀ ' is a N × L|G| matrix and � b g is a L|G| × 1 vector. Thus, the matrix formulation of Eq 2 is:

Elastic net penalty. Because of the large number of predictors as compared to the small number of samples, we use the elastic net penalty, which is a generalization of both ridge and lasso penalties. The elastic net fits the following objective:

Here k � k 1 represents the ℓ 1 -norm and k � k 2 represents the ℓ 2 -norm. For the elastic net, we used the following ranges of hyperparameter values: λ 2 {10 −4 , 10 −3 , . . . , 1}, a 2 {0.1, 0.3, . . ., 0.9}. For lasso, we used λ 2 {10 −5 , . . ., 1}. For ridge, when we used {10 −5 , . . ., 1}, we found that the optimal value selected in some cases was the maximum value of λ = 1. We thus expanded the range to {10 −5 , . . ., 10 6 } to ensure that we were not missing better hyperparameters at larger values. At this point, the optimal λ was found to be 100.

Hyperparameter tuning. Hyperparameters were selected using leave-one-out cross-validation (LOOCV). The hyperparameter (or pair of hyperparameters, for elastic net) that minimizes the mean-squared error on the held-out observations is selected. More specifically, we first fix a hyperparameter (λ, a). Then, for a given gene g and row index i, extract the ith row of X g t and X G tÀ ' . We refer to this extracted validation set as ðX g t Þ i (target) and ðX g tÀ ' Þ i (predictors). The remaining data is the training set, ðX g t Þ À i (target) and ðX G tÀ ' Þ À i (predictors). First, letb g ðl;aÞ;i be theb g ELASTIC NET that is fit from the training set.

b g ðl;aÞ;i ¼ arg min

We then compute prediction error on the validation set, kðX g t Þ i À ðX G tÀ ' Þ ib g ðl;aÞ;i k 2 2 Þ. We repeat the fitb g ðl;aÞ;i and error for every row index i of X g t and for every gene g. The mean heldout cross-validation error for (λ, a) is:

The (λ, a) that minimizes the error in Eq 13 is selected. Permuted coefficients. We evaluate the significance of any given edge g 0 ! g through permutation. In detail, we remove the time dependency between g 0 and g via permutations of individual gene temporal profiles over time.

We first generate a single permuted data setX g t . For each gene, we independently shuffle the temporal profile of each gene g 2 {1, . . ., |G|} across time (Fig 1A) . This is done separately for distinct replicates. We wish to model the hypothesis of no causal relations from any gene g 0 2 ¬g, upon a given effect gene g. We use the unpermuted values of the effect gene X g t and the permuted values of all other causal genes g 0 2 ¬g, asX :g t . The effect gene g remains unpermuted, as we do not consider self-regulatory loops.

Permutation-based causal coefficientsb g 0 ;g ' are then fit as

We use these coefficients to perform FDR calibration. Edge FDR. The result of the elastic net VAR model is a complete network whose edges are weighted according to the estimated regression coefficients.

For each lag ℓ 2 {1, . . ., L} and effect gene g, we control the edge FDR at �0.05 by finding the threshold T g ' such that P g 0 2:g 1fjb g 0 ;g ' j > T g ' g P g 0 2:g 1fjb g 0 ;g ' j > T g ' g þ

For each gene pair (g 0 , g), g 0 2 ¬g, a directed edge g 0 ! g exists if, for at least one of the lags ℓ 2 {1, . . ., L}, jb g 0 ;g ' j > T g ' . Stability selection. Stability selection is used to ensure the robustness of BETS to small sample size. Stability selection is a method for high-dimensional graph estimation that uses bootstrap samples [82] . While the authors prove finite sample control for the family-wise error rate (FWER), we are interested in controlling the false discovery rate (FDR). This procedure draws B = 1000 bootstrap samples, where each sample consists of N rows drawn with replacement from output X g t and input X G tÀ ' (Eq 10). Let the jth bootstrap sample from the original data be X g;j t and X G;j tÀ ' . A set of N row indices, I j , are sampled with replacement from [1, 2, . . ., N]. X g;j t and X G;j tÀ ' are created by choosing the rows I j of X g t and X G tÀ ' . X g;j t is an N × 1 vector and X G;j tÀ ' is an N × L|G| matrix. Now consider the permuted output and input,X g t andX G tÀ ' , constructed fromX g t (Eq 10). Let the jth bootstrap sample from the permuted data beX g;j t andX G;j tÀ ' .X g;j t andX G;j tÀ ' are created by choosing the rows I j ofX g t andX G tÀ ' . Thus, the jth bootstrap sample for both the original and permuted data sets use the same row indices I j .

For each of the 1000 bootstrap samples, we infer a network using the elastic net fit and edge FDR procedure described earlier. Each edge g 0 ! g's selection frequency, π g 0 , g (the frequency of g 0 ! g among the 1000 bootstrap networks) is computed (Fig 1B) .

Stability FDR. To determine the appropriate cutoff for the selection frequency of each edge (π g 0 , g ), we generate a null distribution of selection frequencies using permutations. First, we generate a second permuted data setX g t in which we again independently shuffle the temporal profile of each gene g 2 {1, . . ., |G|} across time. This is done separately for distinct replicates.

We run the stability selection procedure onX g t as if it were X g t , using the same set of row indices I j to generate the bootstrap samples, and usingX g t to generate the permuted coefficients.

After running for all B = 1000 bootstrap samples, we obtain the null selection frequency of each edge,p g 0 ;g .

We control the stability FDR at 0.2 by finding the threshold T b such that P g 0 2:g 1fp g 0 ;g > T b g P g 0 2:g 1fp g 0 ;g > T b g þ P g 0 2:g 1fp g 0 ;g > T b g � 0:2:

Because the maximum lag is 2, each edge g 0 ! g has two possible lags and thus two selection frequencies. The lag with larger absolute value of average coefficient across the 1, 000 networks is considered in both the permuted and the real empirical distributions. So, if jb g 0 ;g 1 j exceeds jb g 0 ;g 2 j, the lag is said to be 1 and the selection frequency p 1 g 0 ;g is used. Network inference performance metrics. Refer to every network edge inferred by a method as a positive and every missing edge as a negative. Let TP be True Positives, FP be False Positives, TN be True Negatives, and FN be False Negatives. Let TPR be True Positive Rate, (i.e., recall), and FPR be False Positive Rate. Then, we have

In the DREAM benchmark, each network inference method is evaluated by comparing the true network (i.e., the network used to generate the synthetic data) with the inferred network at different thresholds for edge inclusion. The two main evaluation metrics are Area Under the Receiver Operating Characteristic curve (AUROC) and Area Under the Precision-Recall curve (AUPR). AUROC plots TPR on the y axis and FPR on the x axis. AUPR plots precision on the y axis and recall on the x axis. When the number of negatives greatly exceeds the number of positives, as with gene networks, which are typically sparse, AUPR is a more relevant metric [83] .

BETS is available for download on Github at https://github.com/lujonathanh/BETS. The software is licensed under the terms of the Apache License, version 2.0. The analysis code is available at https://zenodo.org/record/4009546#.X5XEh0JKg1g.

DREAM Network Inference Challenge. There were five data sets in the DREAM4 Network Inference Challenge, each consisting of ten time series of 21 time points and 100 genes [45, 84] . For the first half of the time series, a "drug perturbation" was applied; this affected about 1/3 of genes. For the second half, the perturbation was removed and the system was allowed to relax back to the wild-type state.

Glucocorticoid gene expression data. We analyzed RNA-sequencing data from a set of experiments developed to study glucocorticoid receptors (GRs) in the human adenocarcinoma and lung model cell line, A549 [6] . There was an original exposure data set of 4 replicates in which cells were stimulated by the glucocorticoid dexamethasone (dex), and gene expression was profiled at {0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12} hours of dex stimulation. There was also an unperturbed data set of 3 replicates in which cells were exposed to dex for 12 hours, after which the conditioned media was replaced and dex removed. Gene expression was profiled at {0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12} hours after dex removal. We integrated the original exposure and unperturbed data into a joint data set with 7 replicates. The original exposure data set is available at the Gene Expression Omnibus (GEO), with reference numbers listed for the rows that list "RNASeq" as the Assay under the column "Experiment_GEO_Series" in Supplementary Table 3 of [6] . The GEO accession numbers for time points {0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12} of the original exposure data set are GSE91305, GSE91198, GSE91311, GSE91358, GSE91303, GSE91243, GSE91281, GSE91229, GSE91255, GSE91284, GSE91222, and GSE91212, respectively. The unperturbed data set is available at GEO accession number GSE144662.

We selected 2768 genes for analysis, which had average expression > 2 Transcripts Per kilobase Million (TPM) and were differentially expressed in the original exposure data. A gene was called differentially expressed if its expression at any time point differed from its expression at time 0, ascertained by running edgeR (FDR � 0.05) [6] . We added NR3C1, which encodes the glucocorticoid receptor (GR). NR3C1 was not found to be differentially expressed at FDR � 0.05.

After genes were selected, gene expression TPM were log-normalized and corrected for surrogate variables using SVAseq [85] . Each gene's temporal profile was centered to have mean zero across time. In the original exposure data, all replicates besides replicate 1 had a measurement for each time point. Replicate 1 was missing time points 5 and 6 hrs, so we imputed these values using a linear interpolation from time points 4 and 7 hrs in the log-transformed, surrogate-corrected space.

Overexpression transcriptional time-series data. There were ten overexpression data sets, in which each of the transcription factors CEBPB, CEBPD, FOSL2, FOXO1, FOXO3, KLF6, KLF9, KLF15, POU5F1, and TFCP2L1 was separately overexpressed across 12 hours of dex stimulation. Each overexpression data set had three replicates; gene expression was profiled after {0, 1, 4, 8, 12} hours of dex stimulation. The same 2768 genes were selected and the same normalization and SVAseq correction as earlier was performed. The overexpression data sets are available at the Gene Expression Omnibus at GEO accession number GSE144660.

Application of methods to the data DREAM benchmarking. We ran the methods BETS, Enet, CSId [44] , Jump3 [36] , CLR [27] , MRNET [28] , ARACNE [29] , SWING-RF [51] , and SWING-Lasso [51] on the DREAM challenge. In BETS, inferred edges were ranked by their selection frequency for calculating AUPR and AUROC. In Enet, edges were ranked by the absolute value of their coefficient. The Python3 version of CSId was run after obtaining it from correspondence with Dr. Penfold.

Jump3 required setting the "systematic noise" and "observational noise" parameters. We used Dr. Huynh-Thu's settings on the DREAM challenge, with systematic noise at 1e − 4 and observational noise at 0.01 times the value of the gene's expression. ARACNE, MRNET, and CLR were run using the minet R library. BETS, Enet, CSId, and Jump3 were run on a single node without parallelization. The node had 28 cores, 128 GB of memory, and 2.4 GHz processor speed. ARACNE, MRNET, and CLR were run on a 4 GB RAM, Intel Core i5 1.3 GHz laptop.

Network analysis: Gene annotations. We considered genes with three possible labels: immune system, metabolism, or transcription factor. Immune genes were labeled as such using two sources. The first source is the Gene Ontology (GO) annotation "Immune" (GO:0002376) [52] . We applied this label when the evidence codes were one of EXP, IDA, IGI, IMP, IPI, IC, TAS. The second source is the Gene Ontology Consortium's curated, ranked list of immune-related genes based on multiple databases and experimental evidence [54] . For the GO annotation, we selected all genes with score � 7. This resulted in 616 immune genes overall, and 109 immune genes in our list of 2768 genes.

Metabolic genes were called using two sources. The first source is the GO annotation "carbohydrate metabolic process" GO:0005975 [52] . We applied this label when the evidence codes were EXP, IDA, IGI, IMP, IPI, IC, TAS. The second source is the Gene Set Enrichment Analysis (GSEA)-curated list of metabolic-related genes [53] . We searched only among those with experimental evidence: the Canonical, KEGG, BIOCARTA, and Reactome pathways. We used the following four search queries: "gluconeogenesis OR (glucose AND metabolism) OR glycolysis," "lipid AND metabolism," "Diabetes," "Obesity." This resulted in 544 metabolic genes overall, of which 120 were in our gene list. 65 genes were both immune and metabolic overall; 12 of these were in our gene list.

Transcription factors (TFs) were called using the Bioguo database of human TFs [55] . There were 1463 TFs overall, of which 226 were present in our gene list.

Experimental interactions. We created a list of experimentally validated interactions from the BIOGRID Homo sapiens Protein-Protein Interactions database [56] . Proteins were mapped to genes using BioMart from Ensembl 94 [86] . Among genes in our gene list, there were 17, 990 BIOGRID interactions.

Validation on overexpression data. The overexpression data had four time points with 1 to 4 hour time gaps, unlike the original 12 time points with 0.5 to 2 hour time gaps. On the overexpression data, we used a VAR model that regressed each effect gene's expression level on its previous expression level and the causal gene's previous expression level, assuming normal noise � t � N ð0; 1Þ:

No regularization was included, and ordinary least squares was used to fit the equation. The expression X g 0 tÀ 1 of a causal gene g 0 is fit as a single predictor without the other expression. Lag 1, not 2, is used due to the larger time gaps.

Validation on lung trans-eQTLs in GTEx v6. Trans-eQTLs were discovered using the Genotype Tissue Expression (GTEx) v6 data [14, 63] . First, we mapped our genes from hg38 to hg19. For every edge g 0 ! g, we tested the set of genetic variants within 20 kilobases of g 0 for trans-eQTL association with g [87] . Specifically, we computed the p-value for linear association of each variant with the corresponding effect gene g using MatrixEQTL [88] . A null distribution was generated by taking every edge g 0 ! g, permuting the effect gene g's expression values, and repeating the linear association test. FDR over test statistics was calculated using q-value [89] . Because not every causal gene g 0 had a cis-eQTL, only 26,839 edges (84% of the original 31,945 edges) were tested. Table. (PDF) S1 Table. DREAM4 100-gene network inference results, AUPR. DBN is dynamic Bayesian network, DT is decision tree, GP is Gaussian process, MI is mutual information, ODE is ordinary differential equation, VAR is vector autoregression. The references that reported ebdbnet, ScanBMA, and LASSO did not provide AUPR values for individual networks. Algorithms that were run in-house were ARACNE, BETS, CLR, CSId, Enet, Jump3, MRNET, SWING-Lasso, and SWING-RF. Where reported literature values were available, they were consistent with these values. Values for CSIc, G1DBN, GCCA, GP4GRN, TSNI, VBSSMa and VBSSMb were taken from [49] . Values for ebdnet, LASSO and ScanBMA, were taken from [40] . Values for dynGENIE3, GENIE3, OKVAR-Boost and tl-CLR were taken from [37] . Values for Inferelator and Jump3 were taken from [36] . Related to Fig 2. (DOCX) S2 Table. DREAM4 100-gene network inference results, AUROC. DBN is dynamic Bayesian network, DT is decision tree, GP is Gaussian process, MI is mutual information, ODE is ordinary differential equation, VAR is vector autoregression. The references that reported ebdbnet, ScanBMA, and LASSO did not provide AUROC values for individual networks. Algorithms that were run in-house were ARACNE, BETS, CLR, CSId, Enet, Jump3, MRNET, SWIN-G-Lasso, and SWING-RF. Values for CSIc, G1DBN, GCCA, GP4GRN, TSNI, VBSSMa and VBSSMb were taken from [49] . Values for ebdnet, LASSO, and ScanBMA, were taken from [40] . Related to Fig 2. (DOCX) S3 Table. Results of in-house algorithms on DREAM4 100-gene network inference. AUPR, AUROC, and Time indicate average AUPR, AUROC, and time over the five networks, respectively. BETS and Enet are in bold to indicate that they are our own developed methods, based on vector autoregression. SWING-RF [51] and Jump3 [36] are decision tree methods. CSId is a Gaussian process method [44] . CLR [27] , MRNET [90] , and ARACNE [29] are mutual information methods. SWING-Lasso is a vector autoregression method [51] . Related to Fig 2. (DOCX) S4 Table. Improvement on DREAM4 100-gene network inference from bootstrap. For each AUROC or AUPR column, the average is the listed value and the standard deviation is listed in parentheses. "Coefficient" denotes the result when ranking edges by their fitted coefficient, as in the original method. "Bootstrap" denotes the results when ranking edges by the frequency by which they appear in the bootstrap networks. 

Studying and modelling dynamic biological processes using timeseries gene expression data

Bayesian factor regression models in the "large p, small n" paradigm

High-dimensional statistics with a view toward applications in biology

Circadian clock function in Arabidopsis thaliana: time beyond transcription

Learning non-stationary dynamic Bayesian networks

Glucocorticoid receptor recruits to enhancers and drives activation by motif-directed binding

Immune regulation by glucocorticoids

World first coronavirus treatment approved for NHS use by government

COVID-19 Treatment Guidelines: Corticosteroids

World Health Organization. Corticosteroids for COVID-19

Infectious Diseases Society of America Guidelines on the Treatment and Management of Patients with COVID-19

Mechanisms of glucocorticoid-induced insulin resistance: focus on adipose tissue function and lipid metabolism

The glucocorticoid contribution to obesity

Genetic effects on gene expression across human tissues

Reverse Engineering of Genome-wide Gene Regulatory Networks from Gene Expression Data

Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process

Grouped graphical Granger modeling for gene expression regulatory networks discovery

Reconstructing Causal Biological Networks through Active Learning

Estimating high-dimensional intervention effects from observational data. The Annals of Statistics

Active Learning of Causal Bayes Net Structure

Joint estimation of causal effects from observational and intervention gene expression data

Two optimal strategies for active learning of causal models from interventional data

Active learning of causal networks with intervention experiments and optimal designs

An introduction to Gaussian Bayesian networks

DREAM3: network inference using dynamic context likelihood of relatedness and the inferelator

Experimental assessment of static and dynamic algorithms for gene regulation inference from time series expression data

Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles

Information-theoretic inference of large transcriptional regulatory networks

ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context

Reverse engineering of gene networks from time-course data by an information theoretic approach

Testing for causality

Inference of gene regulatory networks and compound mode of action from time course gene expression profiles

The Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo

Induction of decision trees

Classification and regression trees. Routledge

Combining tree-based and dynamical systems for the inference of gene regulatory networks

dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data

Inferring dynamic genetic networks with low order independencies

Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks

Fast Bayesian inference for gene regulatory networks using ScanBMA

A Bayesian approach to reconstructing genetic regulatory networks with hidden factors

An empirical Bayesian method for estimating biological networks from temporal microarray data

Learning gene regulatory networks from gene expression measurements using non-parametric molecular kinetics

CSI: a nonparametric Bayesian approach to network inference from multiple perturbed time series gene expression data

Generating realistic in silico gene networks for performance assessment of reverse engineering methods

Sparse inverse covariance estimation with the graphical lasso

Regularization and variable selection via the elastic net

Regression shrinkage and selection via the lasso

How to infer gene networks from expression profiles

Inferring regulatory networks from expression data using treebased methods

Windowed Granger causal inference strategy improves discovery of gene regulatory networks

Gene Ontology: tool for the unification of biology

Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

Gene Ontology Consoritum's Curated List of Immune Genes

AnimalTFDB: a comprehensive animal transcription factor database

The BioGRID interaction database: 2017 update

SOCS-1 and SOCS-3 block insulin signaling by ubiquitin-mediated degradation of IRS1 and IRS2

Suppressor of cytokine signaling-3 provides a novel interface in the cross-talk between angiotensin II and insulin signaling systems

Suppressor of cytokine signaling (SOCS) 1 regulates IL-4-activated insulin receptor substrate (IRS)-2 tyrosine phosphorylation in monocytes and macrophages via the proteasome

A novel pathway for vitamin A signaling mediated by RXR heterodimerization with NGFI-B and NURR1

Orphan receptor TR3 attenuates the p300-induced acetylation of retinoid X receptor-α

Distant regulatory effects of genetic variation in multiple human tissues

PPARγ regulates adipocyte cholesterol metabolism via oxidized LDL receptor 1

Upregulation of OLR1 and IL17A genes and their association with blood glucose and lipid levels in femoropopliteal artery disease. Experimental and Therapeutic Medicine

Oxidized LDL receptor 1 gene polymorphism in patients with metabolic syndrome

LOX-1 boosts immunity

C-type lectin-like receptor LOX-1 promotes dendritic cell-mediated class-switched B cell responses

Inflammation-induced effector CD4+ T cell interstitial migration is alpha-v integrin dependent

Snai2 is a new target to mediate glucocorticoid signaling on breast cancer cell migration

The SHP-1 protein tyrosine phosphatase negatively modulates glucose homeostasis

Deficient SOCS3 and SHP-1 expression in psoriatic T cells

Macrophages of multiple sclerosis patients display deficient SHP-1 expression and enhanced inflammatory phenotype

KEGG: new perspectives on genomes, pathways, diseases and drugs

A continuous tumor-cell line from a human lung carcinoma with properties of type II alveolar epithelial cells

On the dependency of cellular protein levels on mRNA abundance

Advantages and limitations of current network inference methods

Revealing strengths and weaknesses of methods for gene network inference

Utility and limitations of using gene expression data to identify functional associations

Towards inferring causal gene regulatory networks from single cell expression measurements

Stability selection

The relationship between Precision-Recall and ROC curves

The DREAM4 In-silico Network Challenge: Training data, gold standards, and supplementary information

Capturing Heterogeneity in Gene Expression Studies by Surrogate Variable Analysis

Context-specific and differential gene coexpression networks via Bayesian biclustering models

Matrix eQTL: ultra fast eQTL analysis via large matrix operations

Statistical significance for genomewide studies

minet: AR/Bioconductor package for inferring large transcriptional networks using mutual information

The authors would like to thank Gregory Darnell, Derek Aguiar, Ariel Gewirtz, Allison Chaney, Isabella Grabski, Cristina Anastase, and Genna Gliner for helpful discussion, feedback, and generosity in running cluster jobs; and Jian Peng for productive discussion and helpful comments.The authors gratefully acknowledge that this work was performed using the Princeton Research Computing resources sponsored by the Princeton Institute for Computational Science and Engineering (PICSciE) at Princeton University.

Conceptualization: Bianca Dumitrascu, Barbara E. Engelhardt.