Gene regulatory network inference is essential to uncover complex relationships among gene pathways and inform downstream experiments, ultimately paving the way for 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 but additionally infers whether the causal effects are activating or inhibitory. We apply BETS to transcriptional time series data of 2, 768 differentially-expressed genes from A549 cells exposed to glucocorticoids over a period of 12 hours. We identify a network of 2, 768 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 freely 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 2 statistical models of gene regulatory dynamics [1]. These time series data present a unique opportunity to 3 use the coordinated transcriptional response to environmental exposure to infer causal relationships between 4 genes. However, there are several challenges to overcome in the analysis of time series transcriptomic data. 5 These data are generally high-dimensional: the number of quantified gene transcripts-approximately 20,000 6 in human samples-often dramatically exceeds the number of available time points and samples. Many 7 classical statistical assumptions fail to hold in this high-dimensional regime [2, 3]. Moreover, the large 8 number of gene transcripts poses a computational burden, as the number of possible edges in a gene network 9 grows quadratically. Finally, a transcriptional time series often has a small number of time points, and 10 those time points are often not uniformly spaced; furthermore, because transcriptional time series data often 11 quantify transcription post exposure, the time series is not stationary, and genes respond to the exposure 12 and return to baseline at different rates [4, 5]. 13 In this work, we develop an approach that uses the gene transcriptional time series following glucocorticoid 14 (GC) exposure to build a directed gene network...