Inferring the kinetics of stochastic gene expression from singlecell RNAsequencing data
 Jong Kyoung Kim^{1} and
 John C Marioni^{1}Email author
DOI: 10.1186/gb2013141r7
© Kim and Marioni licensee Springer. 2013
Received: 7 November 2012
Accepted: 28 January 2013
Published: 28 January 2013
Abstract
Background
Genetically identical populations of cells grown in the same environmental condition show substantial variability in gene expression profiles. Although singlecell RNAseq provides an opportunity to explore this phenomenon, statistical methods need to be developed to interpret the variability of gene expression counts.
Results
We develop a statistical framework for studying the kinetics of stochastic gene expression from singlecell RNAseq data. By applying our model to a singlecell RNAseq dataset generated by profiling mouse embryonic stem cells, we find that the inferred kinetic parameters are consistent with RNA polymerase II binding and chromatin modifications. Our results suggest that histone modifications affect transcriptional bursting by modulating both burst size and frequency. Furthermore, we show that our model can be used to identify genes with slow promoter kinetics, which are important for probabilistic differentiation of embryonic stem cells.
Conclusions
We conclude that the proposed statistical model provides a flexible and efficient way to investigate the kinetics of transcription.
Keywords
gene regulation RNAseq singlecell statistics transcriptional burstBackground
RNAsequencing (RNAseq) is a recently developed approach that allows an unbiased examination of the transcriptome to be performed using highthroughput DNA sequencing [1–3]. Compared to gene expression microarrays, the previous gold standard for genomewide quantification of gene expression levels, RNAseq has some specific advantages: it allows splicing to be assayed in an unbiased manner [4], it better enables the measurement of expression levels over a wide dynamic range [1], and it allows allelespecific expression to be interrogated [5, 6].
Until recently, most RNAsequencing experiments began with a large population of cells (> 10^{5}), and, as a result, the gene expression counts obtained can be viewed as an average across that population. However, recent developments in sequencing technology have enabled the use of much smaller volumes of starting material, and several groups have described protocols for assaying the transcriptome of single cells [7–11]. This is vital in many biological contexts, such as early embryonic development and tumor etiology, where it is expected that different cells will have distinctive expression profiles. Furthermore, even in tissues that are typically considered to consist of homogeneous populations of cells, intercellular variability in gene expression levels can be considerable. For example, the cells of a genetically identical population grown in the same environment have been shown to display substantial variability in the total number of mRNA molecules that they contain [12–14]. This variability can be partially explained by noting that gene expression levels are regulated by combinatorial interactions between numerous cellular components, where these interactions involve random biochemical reactions [12, 13, 15].
More generally, singlecell imaging methods (e.g., RNA fluorescence in situ hybridization or FISH) have been widely applied to elucidate the principles of gene expression regulation in vivo [16]. These studies have observed that: i) gene expression is heterogeneous [12–14]; ii) genes fluctuate between an 'on' and 'off' promoter state and transcripts are produced in bursts [17–19]; and iii) the transition to the 'on' state requires multiple ratelimiting steps that are determined by many sequential interactions between regulators and chromatin, but the transition to the 'off' state can be determined by a single ratelimiting step [16]. Some examples of the stochastic processes that play a role in the transition to the 'on' state are the recruitment of nucleosome remodelers and histonemodifying enzymes by activators, the rate at which RNA polymerase II (PolII) escapes from the core promoter to produce short RNA molecules prior to pausing, and the rate at which PolII leaves pausing and enters productive elongation [15].
One situation where stochastic fluctuation in gene expression levels plays an important role is in the regulation of mouse embryonic stem (ES) cells [14]. Mouse ES cells are derived from the inner cell mass (ICM) or the epiblast of the preimplantation blastocyst [20], and they can proliferate in the same undifferentiated state indefinitely whilst retaining the ability to differentiate into all adult cell lineages. These two hallmarks of ES cells are conferred by tightly controlled gene regulatory networks [21]. However, growing evidence suggests that the ability of an individual ES cell to differentiate into an adult cell type at a specific time is determined stochastically [14, 22]. In particular, the expression levels of key regulatory genes, such as Nanog, Stella, and Rex1, which are markers of pluripotency, are heterogeneous in ES cells even though the cells are cultured in the same condition [23]. This implies that ES cells exist in a dynamic equilibrium between states that show different propensities for differentiation [22–24].
Here, we develop a statistical framework motivated by a kinetic model for transcriptional bursting to model the biological variability present in singlecell RNAseq data. The framework derived makes it easy to perform parameter fitting and allows the kinetics of transcription to be investigated. We apply our model to singlecell RNAsequencing data generated from mouse ES cells and demonstrate that the estimated parameters are consistent with promoter kinetics inferred from RNA polymerase II binding and chromatin state profiles.
Results
A kinetic model for stochastic gene expression
Moreover, when the gene is in the active promoter state, it is assumed to be transcribed at a rate, s, per unit time and the number of mRNA molecules of the gene is assumed to decay at a rate, d, per unit time. Subsequently, transcriptional bursting can be characterized by two parameters: the average number of synthesized mRNA molecules while a gene remains in an active state (burst size or transcriptional efficiency, s/k_{off}) and the frequency at which bursts occur per unit time (burst frequency, k_{on}) [18, 26–28].
Given these four kinetic parameters, a set of differential equations has been derived describing how the number of mRNA molecules of a given gene within a cell, x, changes over time (Additional file 1; [17]). The steady state distribution of these equations has been shown to take the form [17–19]:
As noted previously, the four kinetic parameters are all measured in units of time. However, since the inverse of the decay rate, 1/d, denotes the average lifetime of an mRNA molecule, it can be used to normalize the other kinetic parameters so that they are independent of time [16, 17]. This is equivalent to setting d = 1 in (2), and we do this henceforth.
A Poissonbeta model
The parameters of the steadystate solution (2) have previously been estimated from observed data using two different approaches. The first approach is to match the first three moments to their empirical values [17]. Although this method is straightforward and computationally efficient, it does not guarantee that the estimates are within the parameter space. To overcome this problem the maximum likelihood estimates of the parameters can be found using a numerical optimization approach [18]. However, the computation of the confluent hypergeometric function, _{1}F_{1}, is difficult, because there is no numerical method for its accurate, fast and reliable computation within all parameter values [29]. Furthermore, when the number of observations is small (less than 100), the maximum likelihood approach sometimes gives unrealistically large estimates of the kinetic parameters [18].
where p is an auxiliary variable following a beta distribution. The marginal distribution P(xs, k_{on}, k_{off}), which is known as the Poissonbeta distribution (PoBe) [30], takes the same form as the steadystate distribution described in equation (2).
Interestingly, the mean of the auxiliary variable p is equal to the fraction of time that a gene spends in the active state (1). Further, Smiley and Proulx [31] showed that if a gene's expression level oscillates between 0 and s/d and the maximum expression level of the gene (s/d) is set to 1, then given the twostate model, its stationary distribution takes the same form as that of p.
Given count measurements from RNAsequencing data, we assume that the number of reads mapped to a gene is proportional to the expression of the relevant mRNA molecule in the cell under study, and thus the parameters of the kinetic model can be inferred using a Bayesian hierarchical approach, such as a Gibbs sampler.
One of the most significant challenges in applying the kinetic model for gene expression is interpreting the parameters. As noted in a recent review by Munsky et al. [25], when k_{on} and k_{off} are large the transitions between the promoter states are rapid, resulting in a Poisson or negative binomiallike distribution of the number of mRNA molecules [18, 19]. In the context of fitting the kinetic model to real data this corresponds to areas of the parameter space where the three parameters are not identifiable (Method; Figure 1B). By contrast, when k_{on} and k_{off} are small, there are relatively few transitions between the two promoter states and the resulting distribution of gene expression molecules between different cells is bimodal  here all three parameters are identifiable (Figure 1BE).
In practice, to ensure that the parameters are statistically identifiable, we suggest fitting three models (Poisson, negative binomial and Poissonbeta) to each gene before using a goodnessoffit statistic to determine whether there is evidence that the parameters of the Poissonbeta model can be identified unambiguously (Methods). An alternative approach would be to fit a hierarchical Bayesian model to each gene and to use this to determine the best fitting distribution.
Assessing the reliability of the Poissonbeta model
Singlecell RNAsequencing was recently used to assay the transcriptome of 12 mouse ES cells derived from the ICM at embryonic day 3.5 (E3.5) [8]. To explore the transcriptional kinetics of ES cells, we fitted the Poissonbeta model to these data (Methods). Before interpreting the inferred kinetic parameters, it is necessary to: i) account for the high amount of technical variability present in singlecell RNAseq data; ii) consider whether the parameter estimates are statistically identifiable; and iii) assess whether we can draw meaningful inferences about transcriptional kinetics based on gene expression measurements from 12 cells.
Accurately quantifying the technical variability present in singlecell RNAseq data is challenging. While experimental approaches vary, most suggest that when replicate libraries are generated from small quantities of RNA (taken from the same, large, population of RNA), the resulting read counts display more technical variability, especially for lowly expressed genes, than is observed in populationbased RNAsequencing analyses [7, 10, 11]. This is likely due to experimental factors such as the efficiency of the RT step and the PCR amplification when small quantities of starting material are considered [7, 10, 11]. Some attempts have been made to characterize the technical variation using spikeins [11] but evidence for the efficacy of such approaches is still limited. Given these challenges and the limitations of current experimental approaches, we instead removed lowly expressed genes that are most likely to display high technical variability [7, 10, 11]. We considered a gene as lowly expressed if the maximum normalized read count was less than 50. This cutoff was determined using technical replicate data generated using the same protocol applied to the 12 ES cells [8] or oocytes [7] (Additional file 1, Figure S2). Across the set of 18,735 genes that were expressed in at least one cell, 12,551 genes had an expression level above this cutoff, and we fitted the Poissonbeta model separately to each of these genes (Additional file 1, Figure S3).
 1.
Genes with relatively large values of k_{off,i}and low values of k_{on,i}. This corresponds to genes that have a low expression count in most cells and high expression in a small number of cells (typically one). When we simulated data from the Poissonbeta model with parameter values in this range (Methods), we found that k_{on,i}was estimated accurately, but that both k_{off,i}and s_{ i } were underestimated (Additional file 1, Figures S6, S11).
 2.
Genes with large values of k_{off,i}and k_{on,i}(Additional file 1, Figure S3). This set of genes are typically highly expressed (Additional file 1, Figure S12C, F) with a relatively low amount of variability across cells, as evidenced by the large values of k_{off,i}and k_{on,i}. While it is possible that this set of genes do have very fast promoter kinetics, statistically it is impossible to distinguish this from their being permanently in an active state (that is, k_{off,i}is equal to or very close to zero). Hence, it is impossible to interpret either the raw or the derived parameters in this situation (Additional file 1, Figure S12). More generally, this explains why we do not observe many low values of k_{off,i}in our set of parameter estimates. Moreover, it helps to explain some of the identifiability problems that other approaches have encountered when estimating k_{off,i}.
Given these observations, we focus henceforth on the 10,298 genes that have identifiable estimates of the kinetic parameters. However, before going on to make biological inferences based upon these parameters it is first necessary to assess whether any meaningful conclusions can be drawn from fitting the Poissonbeta model to data from only 12 independent ES cells.
To do this, we fitted the Poissonbeta model to data simulated using the estimated parameters by increasing the number of cells from 3 to 100 (Methods). As expected, the correlation between the parameter estimates and the true values improved as the number of cells increased (Additional file 1, Figure S4S8), with a good agreement when 12 cells were considered (Additional file 1, Figure S6S8). Our simulations also displayed a tendency to underestimate s_{ i }; the extent of the underestimation decreased as the number of cells increased (Additional file 1, Figure S8). One effect of this is a slight bias in the estimated values of k_{on,i}and k_{off,i}(Additional file 1, Figures S6S8). This is not unexpected since s_{ i } can be considered to represent the 'maximum' rate of transcription and, especially when the number of cells is small, a cell where a gene is expressed at this 'maximal' value will not be simulated. Nevertheless, our simulations do provide confidence in the fit of the Poissonbeta model when a moderate number of cells (greater than or equal to 12) are considered. However, it is important to acknowledge that drawing strong biological inferences about the kinetic parameters of individual genes from only 12 cells is difficult  hence, in what follows we consider properties of sets of genes with specific values of the kinetic parameters. This will also help mitigate any effect that technical noise in the measurement of gene expression levels will have upon our interpretation of the data.
Transcriptional kinetics of mouse ES cells
To explore how the kinetic parameters provide information about the regulation of transcription we utilized independent information collected from v6.5 mouse ES cells (derived from the ICM at E3.5) on RNA polymerase II (PolII) occupancy and various histone modifications [32–34]. At the global level, the rates of transcription and gene activation are strongly correlated with the average expression level while the rate of gene inactivation displays a more modest correlation (Additional file 1, Figure S13).
Discussion
The Poissonbeta model provides a convenient statistical framework for modeling singlecell RNAseq data and for studying the kinetics of stochastic gene expression. Since the kinetic parameters of individual genes inferred from the small number of cells are likely to be noisy and may be influenced by technical variability, we focused on the summary properties of genes. Importantly, we confirmed that the kinetic parameters derived from the Poissonbeta model are consistent with PolII binding and chromatin modifications using singlecell RNAseq data generated from mouse ES cells. Our results suggest that the chromatin state of genes, defined by H3K4me3, H3K27me3 and H3K36me3 modifications, affects transcriptional bursting by modulating both burst size and frequency, consistent with a recent study that suggested chromosomal location affected these kinetic characteristics [40].
However, while our model has clear advantages, it also has a number of limitations. First, in this manuscript, we do not address the modeling of the technical variability directly, primarily because our understanding of how different experimental characteristics (RT efficiency, PCR amplification, etc.) might contribute to the noise is very limited. Instead, we focus only on genes that are moderately to highly expressed, since previous singlecell RNAseq studies have shown that such genes display less technical variability. However, as our understanding of the technical variability inherent to singlecell RNAseq increases, it will be important to adapt the model presented herein.
Second, in common with most other biochemical models of gene expression, we assume that the rate of transitions to the 'on' state is governed by a single ratelimiting step. While this assumption facilitates the derivation of a closedform solution for the master equations, and thus the implementation of the Poissonbeta model described in this paper, in higher eukaryotes activation requires many sequential steps [15]. However, the limited experimental data about the relative contribution of the different steps justifies the simplified model presented herein.
Third, the three kinetic parameters are currently measured in units of 'per mRNA average lifetime' since they are normalized by the decay rate. To estimate them in units of 'per second', we should directly measure the decay rates of all genes. This can be done by metabolic labeling of RNA with 4thiouridine coupled with massively parallel sequencing [41]. Another improvement would be to measure the number of mRNA molecules directly rather than using the number of reads as a surrogate, which can be done by accurate digital quantification of transcriptome via digital RNAseq [42].
Finally, our model assumes that the transition times and kinetic parameters are identical for the two alleles of each gene. A recent study established that 39% to 51% of heterozygous loci show allelespecific expression when expression patterns are measured in single cells of a twocell embryo [43]. This suggests that the kinetic parameters and transition times of the underlying Markov chain might differ significantly between the two alleles of a gene. Further, Miyanari et al. [44] showed that Nanog is largely expressed from a single allele in ES cells and can transition between alleles randomly. Such variability can be incorporated into our model by measuring the expression of each allele independently (for example by using MMSEQ [45]), and using these measures as the input to the model. However, the mouse ES cell data we analyzed were generated from an inbred population of mice (C57BL/6J) and, as a result, we could not apply this approach. Examining allelespecific variation using the Poissonbeta model provides an interesting avenue for future research.
Conclusions
To summarize, as the singlecell field progresses towards analyzing the transcriptome of large numbers of individual cells in parallel, it will become increasingly important to develop statistical methods that accurately model stochastic gene expression. In this context, we anticipate that the Poissonbeta model presented here, and other similar approaches, will be vital in maximizing the amount of biological insight that can be obtained from these data.
Materials and methods
Properties of the steadystate solution of the master equation
A hierarchical Bayesian model
 1.Draw s_{ i } for gene i from a gamma distribution:${s}_{i}~\mathsf{\text{Gamma}}\left({s}_{i}{\alpha}_{{s}_{i}},{\beta}_{{s}_{i}}\right)$
 2.Draw k_{on,i}for gene i from a gamma distribution:${{k}_{\mathsf{\text{on,}}}}_{i}~\mathsf{\text{Gamma}}\left({{k}_{\mathsf{\text{on,}}}}_{i}{\alpha}_{{k}_{\mathsf{\text{on}}}{,}_{i}},{\beta}_{{k}_{\mathsf{\text{on}}}{,}_{i}}\right)$
 3.Draw k_{off,i}for gene i from a gamma distribution${{k}_{\mathsf{\text{off,}}}}_{i}~\mathsf{\text{Gamma}}\left({{k}_{\mathsf{\text{off,}}}}_{i}{\alpha}_{{k}_{\mathsf{\text{off}}}{,}_{i}},{\beta}_{{k}_{\mathsf{\text{off}}}{,}_{i}}\right)$
 4.Draw p_{ ij } for gene i and cell j from a beta distribution${p}_{ij}~\mathsf{\text{Beta}}\left({p}_{ij}{k}_{\mathsf{\text{on}},i},{k}_{\mathsf{\text{off}},i}\right)$
 5.Draw x_{ ij } from a Poisson distribution${x}_{ij}~\mathsf{\text{Poisson}}\left({x}_{ij}\mathsf{\text{}}{\mathsf{\text{t}}}_{i}{t}_{j}{s}_{i}{p}_{ij}\right)$
where t_{ i } is the length of gene i (the length of the transcripts measured in bp) and t_{ j } is the normalization factor for cell j. We used the scale normalization method of [46] to estimate the normalization factor for each cell. We make an implicit assumption that the number of reads is proportional to the number of mRNA molecules present in a cell.
The graphical model representing this generative process is shown in Additional file 1, Figure S1.
Learning by collapsed Gibbs sampling
Let $\mathcal{X}=\left\{{x}_{ij}\right\}$ be a set of observed read counts, and $\mathcal{P}=\left\{{p}_{ij}\right\}$ be a set of p_{ ij }. We treat the toplevel variables, shown in the graphical model in Additional file 1, Figure S1, $\mathrm{\Psi}=\left\{{\alpha}_{{k}_{\mathsf{\text{on}},i}},{\beta}_{{k}_{\mathsf{\text{on}},i}},{\alpha}_{{k}_{\mathsf{\text{off}},i}},{\beta}_{{k}_{\mathsf{\text{off}},i}},{\alpha}_{{s}_{i}},{\beta}_{{s}_{i}}\right\}$, as fixed hyperparameters. We derive a collapsed Gibbs sampler to infer all unknown variables $\mathrm{\Theta}=\left\{\mathcal{P},\left\{{k}_{\mathsf{\text{on}},i}\right\},\left\{{k}_{\mathsf{\text{off}},i}\right\},\left\{{s}_{i}\right\}\right\}$ given $\mathcal{X}$. In the following, a subscript with a minus sign that is attached to a set of variables means that the variables indexed by the subscript are excluded from the set.
 1.Sampling p_{ ij }$P\left({p}_{ij}{\mathcal{P}}_{ij},\mathcal{X},\mathrm{\Theta}\backslash \mathcal{P},\mathrm{\Psi}\right)\propto P\left({p}_{ij}{k}_{\mathsf{\text{on}},i},{k}_{\mathsf{\text{off}},i}\right)P\left({x}_{ij}{s}_{i},{p}_{ij}\right)$
 2.Sampling k_{on,i}$P\left({k}_{\mathsf{\text{on}},i}\left\{{k}_{\mathsf{\text{on}},i}\right\},\mathcal{X},\mathrm{\Theta}\backslash \left\{{k}_{\mathsf{\text{on}},i}\right\},\mathrm{\Psi}\right)\propto P\left({k}_{\mathsf{\text{on}},i}{\alpha}_{{k}_{\mathsf{\text{on}},i}},{\beta}_{{k}_{\mathsf{\text{on}},i}}\right){\displaystyle \prod _{j=1}^{J}}P\left({p}_{ij}{k}_{\mathsf{\text{on}},i},{k}_{\mathsf{\text{off}},i}\right)$
 3.Sampling k_{off,i}$P\left({k}_{\mathsf{\text{off}},i}\left\{{k}_{\mathsf{\text{off}},i}\right\},\mathcal{X},\mathrm{\Theta}\backslash \left\{{k}_{\mathsf{\text{off}},i}\right\},\mathrm{\Psi}\right)\propto P\left({k}_{\mathsf{\text{off}},i}{\alpha}_{{k}_{\mathsf{\text{off}},i}},{\beta}_{{k}_{o\mathsf{\text{off,}}i}}\right){\displaystyle \prod _{j=1}^{J}}P\left({p}_{ij}{k}_{\mathsf{\text{on}},i},{k}_{\mathsf{\text{off}},i}\right)$
 4.Sampling s_{ i }$P\left({s}_{i}\left\{{s}_{i}\right\},\mathcal{X},\mathrm{\Theta}\backslash \left\{{s}_{i}\right\},\mathrm{\Psi}\right)\propto P\left({s}_{i}{\alpha}_{{s}_{i}},{\beta}_{{s}_{i}}\right){\displaystyle \prod _{j=1}^{J}}P\left({x}_{ij}{s}_{i},{p}_{ij}\right)$
For the hyperparameters, we used the following settings: ${\alpha}_{{s}_{i}}=1$, ${\beta}_{{s}_{i}}={max}_{j}{x}_{ij}$, ${\alpha}_{{k}_{\mathsf{\text{on}},i}}=1$, ${\beta}_{{k}_{\mathsf{\text{on}},i}}=100$, ${\alpha}_{{k}_{\mathsf{\text{off}},i}}=1$ and ${\beta}_{{k}_{\mathsf{\text{off}},i}}=100.$ We chose the empirical Bayes prior on s_{ i } so that it becomes almost uniform across all realistic ranges for the parameter. The priors on k_{on,i}and k_{off,i}were chosen to place substantial probability across the identifiable parameter space. When we used different priors on k_{on,i}and k_{off,i}(${\alpha}_{{k}_{\mathsf{\text{on}},i}}=1$, ${\beta}_{{k}_{\mathsf{\text{on}},i}}=10,000$, ${\alpha}_{{k}_{\mathsf{\text{off}},i}}=1$, ${\beta}_{{k}_{\mathsf{\text{off}},i}}=10,000$ for more diffuse priors and ${\alpha}_{{k}_{\mathsf{\text{on}},i}}=1$, ${\beta}_{{k}_{\mathsf{\text{on}},i}}=10$, ${\alpha}_{{k}_{\mathsf{\text{off}},i}}=1$, ${\beta}_{{k}_{\mathsf{\text{off}},i}}=10$ for more concentrated priors), the inferred kinetic parameters remained similar, except for the large values of the parameters that were penalized by the concentrated prior. These results suggest that our model is relatively insensitive to the choice of priors (Additional file 1, Figure S9, S10).
Bootstrapbased goodnessoffit test
 1.
Given the maximum likelihood estimates ${\theta}_{n}^{\mathsf{\text{dist}}}$, generate n bootstrap samples ${X}_{1,k}^{*},\dots ,{X}_{n,k}^{*}$ from ${F}_{{\theta}_{n}^{\mathsf{\text{dist}}}}\left(\mathsf{\text{dist}}\in \left\{\mathsf{\text{Poisson}},\mathsf{\text{NB}}\right\}\right)$.
 2.
Compute the maximum likelihood estimates ${\theta}_{n,k}^{*,\mathsf{\text{dist}}}$ from the bootstrap samples.
 3.Estimate the empirical distribution function of the bootstrap samples${F}_{n,k}^{*}\left(x\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}}1\left({X}_{i,k}^{*}\le x\right).$
 4.Compute the KS test statistic ${\mathsf{\text{KS}}}_{n,k}^{*,\mathsf{\text{dist}}}$ such that${\mathsf{\text{KS}}}_{n,k}^{*,\mathsf{\text{dist}}}={max}_{i}{F}_{n,k}^{*}\left({X}_{i,k}^{*}\right){F}_{{\theta}_{n,k}^{*,\mathsf{\text{dist}}}}\left({X}_{i,k}^{*}\right).$
For this study, we set n = 1,000 and B = 1,000.
Estimating the kinetic parameters from synthetic data
Given the posterior means of s_{ i }, k_{on,i}and k_{off,i}for each gene, we generated 3, 6, 12, 20 or 100 independent samples from the Poissonbeta distribution. We run the Gibbs sampling algorithm by setting the total number of Gibbs iterations to 10,000, and computed the posterior means by discarding the first half of the samples in each chain as a burnin period.
Gene activity and pause index for GROSeq
To quantify PolII activity at the promoters and gene body regions, we defined three measures based on the number of mapped reads of GROSeq data [49]. First, gene body activity is defined as N/L where N is the number of GROSeq reads mapped from +1 kb of the transcription start site to the end of a gene, and L is the length of the region. Second, gene promoter activity is defined as the maximum count of reads in a 50 bp window, where we took the maximum among all the windows within the ± 1 kb region of transcription start sites. Finally, we defined the pause index as the ratio of the gene promoter activity (divided by 50) to the gene body activity.
Gene ontology analysis using DAVID
To examine whether particular classes of GO biological processes (GOTERM_BP_FAT) are enriched in the top or bottom 3,000 genes sorted by k_{off,i}/ k_{on,i}, we used the DAVID functional annotation clustering tool (the classification stringency was set to 'medium') [36]. By setting up the 10,298 genes as a background, we chose a representative GO term from each annotation cluster with the Benjaminicorrected P value less than 0.05, providing 18 GO terms in total. The results are in Additional file 2, Table S1 (top 3,000) and Additional file 3, Table S2 (bottom 3,000).
Code availability
The MATLAB source code and a compiled version of the same are available in Additional file 4.
Abbreviations
 bp:

base pair
 ES cell:

embryonic stem cell
 GO:

gene ontology
 ICM:

inner cell mass
 PCR:

polymerase chain reaction
 PoBe:

Poissonbeta
 RNA FISH:

RNA fluorescence in situ hybridization
 RNA PolII:

RNA polymerase II
 RT:

reverse transcriptase.
Declarations
Acknowledgements
We thank Simon Anders, Wolfgang Huber, Anestis Touloumis, Nuno Fonseca and all members of the Marioni group for helpful comments.
Authors’ Affiliations
References
 Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNASeq. Nature Methods. 2008, 5: 621628. 10.1038/nmeth.1226.PubMedView ArticleGoogle Scholar
 Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y: RNAseq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Research. 2008, 18: 15091517. 10.1101/gr.079558.108.PubMedPubMed CentralView ArticleGoogle Scholar
 Wang Z, Gerstein M, Snyder M: RNASeq: a revolutionary tool for transcriptomics. Nature Reviews Genetics. 2009, 10: 5763. 10.1038/nrg2484.PubMedPubMed CentralView ArticleGoogle Scholar
 Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, Kingsmore SF, Schroth GP, Burge CB: Alternative isoform regulation in human tissue transcriptomes. Nature. 2008, 456: 470476. 10.1038/nature07509.PubMedPubMed CentralView ArticleGoogle Scholar
 Degner JF, Marioni JC, Pai AA, Pickrell JK, Nkadori E, Gilad Y, Pritchard JK: Effect of readmapping biases on detecting allelespecific expression from RNAsequencing data. Bioinformatics. 2009, 25: 32073212. 10.1093/bioinformatics/btp579.PubMedPubMed CentralView ArticleGoogle Scholar
 Zhang K, Li JB, Gao Y, Egli D, Xie B, Deng J, Li Z, Lee JH, Aach J, Leproust EM, Eggan K, Church GM: Digital RNA allelotyping reveals tissuespecific and allelespecific gene expression in human. Nature Methods. 2009, 6: 613618. 10.1038/nmeth.1357.PubMedPubMed CentralView ArticleGoogle Scholar
 Tang F, Barbacioru C, Wang Y, Nordman E, Lee C, Xu N, Wang X, Bodeau J, Tuch BB, Siddiqui A, Lao K, Surani MA: mRNASeq wholetranscriptome analysis of a single cell. Nature Methods. 2009, 6: 377382. 10.1038/nmeth.1315.PubMedView ArticleGoogle Scholar
 Tang F, Barbacioru C, Bao S, Lee C, Nordman E, Wang X, Lao K, Surani MA: Tracing the derivation of embryonic stem cells from the inner cell mass by singlecell RNASeq analysis. Cell Stem Cell. 2010, 6: 468478. 10.1016/j.stem.2010.03.015.PubMedPubMed CentralView ArticleGoogle Scholar
 Tang F, Lao K, Surani MA: Development and applications of singlecell transcriptome analysis. Nature Methods. 2011, 8: S6S11. 10.1038/nchembio.740.PubMedPubMed CentralGoogle Scholar
 Islam S, Kjällquist U, Moliner A, Zajac P, Fan JB, Lonnerberg P, Linnarsson S: Characterization of the singlecell transcriptional landscape by highly multiplex RNAseq. Genome Research. 2011, 21: 11601167. 10.1101/gr.110882.110.PubMedPubMed CentralView ArticleGoogle Scholar
 Ramskold D, Luo S, Wang YC, Li R, Deng Q, Faridani OR, Daniels GA, Khrebtukova I, Loring JF, Laurent LC, Schroth GP, Sandberg R: Fulllength mRNASeq from singlecell levels of RNA and individual circulating tumor cells. Nature Biotechnology. 2012, 30: 777782. 10.1038/nbt.2282.PubMedPubMed CentralView ArticleGoogle Scholar
 Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theories to phenotypes. Nature Reviews Genetics. 2005, 6: 451464. 10.1038/nrg1615.PubMedView ArticleGoogle Scholar
 Raj A, van Oudenaarden A: Nature, nurture, or chances: stochastic gene expression and its consequences. Cell. 2008, 135: 216226. 10.1016/j.cell.2008.09.050.PubMedPubMed CentralView ArticleGoogle Scholar
 Eldar A, Elowitz MB: Functional roles for noise in genetic circuits. Nature. 2010, 467: 167173. 10.1038/nature09326.PubMedPubMed CentralView ArticleGoogle Scholar
 Fuda NJ, Ardehali MB, Lis JT: Defining mechanisms that regulate RNA polymerase II transcription in vivo. Nature. 2010, 461: 186192.View ArticleGoogle Scholar
 Larson DR: What do expression dynamics tell us about the mechanism of transcription?. Current Opinion in Genetics & Development. 2011, 21: 591599. 10.1016/j.gde.2011.07.010.View ArticleGoogle Scholar
 Peccoud J, Ycart B: Markovian modelling of gene product synthesis. Theoretical Population Biology. 1995, 48: 222234. 10.1006/tpbi.1995.1027.View ArticleGoogle Scholar
 Raj A, Peskin CS, Tranchin D, Vargas DY, Tyagi S: Stochastic mRNA synthesis in mammalian cells. PLoS Biology. 2006, 4: e30910.1371/journal.pbio.0040309.PubMedPubMed CentralView ArticleGoogle Scholar
 Shahrezaei V, Swain PS: Analytical distributions for stochastic gene expression. Proceedings of the National Academy of Sciences, USA. 2008, 105: 1725617261. 10.1073/pnas.0803850105.View ArticleGoogle Scholar
 Young RA: Control of the embryonic stem cell state. Cell. 2011, 144: 940954. 10.1016/j.cell.2011.01.032.PubMedPubMed CentralView ArticleGoogle Scholar
 Huang S: Cell lineage determination in state space: a systems view brings flexibility to dogmatic canonical rules. PLoS Biology. 2010, 8: e100038010.1371/journal.pbio.1000380.PubMedPubMed CentralView ArticleGoogle Scholar
 Martinez AA, Brickman JM: Gene expression heterogeneities in embryonic stem cell populations: origin and function. Current Opinion in Cell Biology. 2011, 23: 17. 10.1016/j.ceb.2010.12.003.View ArticleGoogle Scholar
 Silva J, Smith A: Capturing Pluripotency. Cell. 2008, 132: 532536. 10.1016/j.cell.2008.02.006.PubMedPubMed CentralView ArticleGoogle Scholar
 Canham MA, Sharov AA, Ko MS, Brickman JM: Functional heterogeneity of embryonic stem cells revealed through translational amplification of an early endodermal transcript. PLoS Biology. 2010, 8: e100037910.1371/journal.pbio.1000379.PubMedPubMed CentralView ArticleGoogle Scholar
 Munsky B, Neuert G, van Oudenaarden A: Using gene expression noise to understand gene regulation. Science. 2012, 336: 183187. 10.1126/science.1216379.PubMedPubMed CentralView ArticleGoogle Scholar
 Skupsky R, Burnett JC, Foley JE, Schaffer DV, Arkin AP: HIV promoter integration site primarily modulates transcriptional burst size rather than frequency. PLoS Computational Biology. 2010, 6: e100095210.1371/journal.pcbi.1000952.PubMedPubMed CentralView ArticleGoogle Scholar
 Batenchuk C, StPierre S, Tepliakova L, Adiga S, Szuto A, Kabbani N, Bell JC, Baetz K, Kaern M: Chromosomal position effects are linked to sir2mediated variation in transcriptional burst size. Biophysical Journal. 2011, 100: L56L58. 10.1016/j.bpj.2011.04.021.PubMedPubMed CentralView ArticleGoogle Scholar
 MillerJensen K, Dey SS, Schaffer DV, Arkin AP: Varying virulence: epigenetic control of expression noise and disease processes. Trends in Biotechnology. 2011, 29: 517525. 10.1016/j.tibtech.2011.05.004.PubMedView ArticleGoogle Scholar
 Muller KE: Computing the confluent hypergeometric function, M(a, b, x). Numerishe Mathematik. 2001, 90: 179196. 10.1007/s002110100285.View ArticleGoogle Scholar
 Johnson NL, Kemp AW, Kotz S: Univariate discrete distributions. 2005, WileyView ArticleGoogle Scholar
 Smiley MW, Proulx SR: Gene expression dynamics in randomly varying environments. Journal of Mathematical Biology. 2010, 61: 231251. 10.1007/s002850090298z.PubMedView ArticleGoogle Scholar
 Mikkelsen TS, Ku M, Jaffe DB, Issac B, Lieberman E, Giannoukos G, Alvarez P, Brockman W, Kim TK, Koche RP, Lee W, Mendenhall E, O'Donovan A, Presser A, Russ C, Xie X, Meissner A, Wernig M, Jaenisch R, Nusbaum C, Lander ES, Bernstein BE: Genomewide maps of chromatin state in pluripotent and lineagecommitted cells. Nature. 2007, 448: 553560. 10.1038/nature06008.PubMedPubMed CentralView ArticleGoogle Scholar
 Ku M, Koche RP, Rheinbay E, Mendenhall EM, Endoh M, Mikkelsen TS, Presser A, Nusbaum C, Xie X, Chi AS, Adli M, Kasif S, Ptaszek LM, Cowan CA, Lander ES, Koseki H, Bernstein BE: Genomewide analysis of PRC1 and PRC2 occupancy identifies two classes of bivalent domains. PLoS Genetics. 2008, 4: e100024210.1371/journal.pgen.1000242.PubMedPubMed CentralView ArticleGoogle Scholar
 Min IM, Waterfall JJ, Core LJ, Munroe RJ, Schimenti J, Lis JT: Regulating RNA polymerase pausing and transcription elongation in embryonic stem cells. Genes & Development. 2011, 25: 742754. 10.1101/gad.2005511.View ArticleGoogle Scholar
 Suganuma T, Workman JL: Signals and combinatorial functions of histone modifications. Annual Review of Biochemistry. 2011, 80: 473499. 10.1146/annurevbiochem061809175347.PubMedView ArticleGoogle Scholar
 Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protocols. 2009, 4: 4457.View ArticleGoogle Scholar
 Gu B, Zhang J, Wang W, Mo L, Zhou Y, Chen L, Liu Y, Zhang M: Global expression of cell surface proteins in embryonic stem cells. PLoS ONE. 2010, 5: e1579510.1371/journal.pone.0015795.PubMedPubMed CentralView ArticleGoogle Scholar
 Cui L, Johkura K, Yue F, Ogiwara N, Okouchi Y, Asanuma K, Sasaki K: Spatial distribution and initial changes of SSEA1 and other cell adhesionrelated molecules on mouse embryonic stem cells before and during differentiation. Journal of Histochemistry & Cytochemistry. 2004, 52: 14471457. 10.1369/jhc.3A6241.2004.View ArticleGoogle Scholar
 HemmatiBrivanlou A, Melton D: Vertebrate embryonic cells will become nerve cells unless told otherwise. Cell. 1997, 88: 1317. 10.1016/S00928674(00)81853X.PubMedView ArticleGoogle Scholar
 Dar RD, Razooky BS, Singh A, Trimeloni TV, McCollum JM, Cox CD, Simpson ML, Weinberger LS: Transcriptional burst frequency and burst size are equally modulated across the human genome. Proceedings of the National Academy of Sciences, USA. 2012, 109: 1745417459. 10.1073/pnas.1213530109.View ArticleGoogle Scholar
 Rabani M, Levin JZ, Fan L, Adiconis X, Raychowdhury R, Garber M, Gnirke A, Nusbaum C, Hacohen N, Friedman N, Amit I, Regev A: Metabolic labeling of RNA uncovers principles of RNA production and degradation dynamics in mammalian cells. Nature Biotechnology. 2011, 29: 436442. 10.1038/nbt.1861.PubMedPubMed CentralView ArticleGoogle Scholar
 Shiroguchi K, Jia TZ, Sims PA, Xie XS: Digital RNA sequencing minimizes sequencedependent bias and amplification noise with optimized singlemolecule barcodes. Proceedings of the National Academy of Sciences, USA. 2012, 109: 13471352. 10.1073/pnas.1118018109.View ArticleGoogle Scholar
 Tang F, Barbacioru C, Nordman E, Bao S, Lee C, Wang X, Tuch BB, Heard E, Lao K, Surani MA: Deterministic and stochastic allele specific gene expression in single mouse blastomeres. PLoS ONE. 2011, 6: e2120810.1371/journal.pone.0021208.PubMedPubMed CentralView ArticleGoogle Scholar
 Miyanari Y, TorresPadilla M: Control of groundstate pluripotency by allelic regulation of Nanog. Nature. 2012, 483: 470473. 10.1038/nature10807.PubMedView ArticleGoogle Scholar
 Turro E, Su SY, Goncalves A, Coin LJ, Richardson S, Lewin A: Haplotype and isoform specific expression estimation using multimapping RNAseq reads. Genome Biology. 2011, 12: R1310.1186/gb2011122r13.PubMedPubMed CentralView ArticleGoogle Scholar
 Anders S, Huber W: Differential expression analysis for sequence count data. Genome Biology. 2010, 11: R10610.1186/gb20101110r106.PubMedPubMed CentralView ArticleGoogle Scholar
 Neal RM: Slice sampling. The Annals of Statistics. 2003, 31: 705767. 10.1214/aos/1056562461.View ArticleGoogle Scholar
 Best DJ, Rayner JCW: Goodness of fit for the Poisson distribution. Statistics & Probability Letters. 1999, 44: 259265. 10.1016/S01677152(99)000176.View ArticleGoogle Scholar
 Core LJ, Waterfall JJ, Lis JT: Nascent RNA sequencing reveals widespread pausing and divergent initiation at human promoters. Science. 2008, 322: 18451848. 10.1126/science.1162228.PubMedPubMed CentralView ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.