 Research
 Open Access
 Published:
Inferring the kinetics of stochastic gene expression from singlecell RNAsequencing data
Genome Biologyvolume 14, Article number: R7 (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.
Background
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
The standard kinetic model for gene expression assumes that a gene can fluctuate randomly between 'on' and 'off' promoter states, where mRNA can be transcribed only in the 'on' state [16, 25] (Figure 1A). If a single ratelimiting step determines the rates of transcription and transitions between the two promoter states [16, 17], the fluctuations between the 'on' and 'off' promoter states can be described by a twostate Markov process where k_{on} is the rate (per unit time) at which a gene becomes active and k_{off} is the rate (per unit time) at which the gene becomes inactive. Consequently, 1/k_{off} and 1/k_{on} describe the average waiting time of a gene in the active and inactive states, respectively, and the (average) fraction of time that a gene spends in the active state is:
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].
To overcome these limitations, we propose an auxiliary variable approach. Specifically, we let:
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).
Using the identifiability criteria outlined in the previous section, we determined that 10,298 (82%) of the 12,551 genes had identifiable parameters at a P value threshold of 0.1 (Methods). The genes with nonidentifiable parameters could be split into two broad categories (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).
As expected, we observed that PolII occupancy was positively correlated with the average expression level (Additional file 1, Figure S14) and burst frequency and, less strongly, with burst size (Figure 2). This is true irrespective of whether PolII levels are calculated in the gene body (P < 10^{16} by the Spearman rank test) or in the promoter region (P < 10^{16}) although the correlation was noticeably higher in the gene body comparison (Figure 2, S14 in Additional file 1). However, when we examined the relationship between burst size and the pause index, defined as the ratio of PolII occupancy in the promoter compared to the gene body, we observed no correlation (P = 0.5558; Figure 2G). Moreover, although we found that both burst frequency and the average proportion of time a gene is transcriptionally inactive were significantly correlated with the pause index (P < 10^{16} for k_{on,i}; P < 10^{13} for k_{off,i}/ (k_{off,i}+ k_{on,i}); Figure 2HI), the correlation is low in both cases, providing only very weak evidence that PolII pausing is associated with burst frequency. A more stringent cutoff that filtered out the lowly expressed genes did not change the results (Additional file 1, Figure S15).
Histone modifications can alter chromatin structure, thereby affecting the regulation of gene expression levels [35]. Two of the most widely studied modifications, H3K4me3, which is associated with active promoters, and H3K27me3, which is associated with genes that have repressed expression levels, have previously been positively and negatively correlated with gene expression levels [32]. Our estimated kinetic parameters are consistent with these observations (see figure legend for statistical details, Figure 3AI, S16 in Additional file 1): genes with a H3K4me3 modification have significantly higher rates of transcriptional bursting and frequency than genes with no modification or those with a H3K27me3 modification. The third histone mark, H3K36me3, which is linked to transcriptional elongation and is enriched over the gene body region [32], was strongly associated with both burst frequency and the fraction of time that a gene spends in the inactive state but more weakly with burst size (Figure 3JL).
Finally, we used gene ontology (GO) analysis [36] to interrogate the set of genes that showed characteristics associated with transcriptional bursting (that is, a low value of k_{on,i}and a relatively high value of k_{off,i}). To do this, we sorted all 10,298 genes in descending order according to the ratio of k_{off,i}to k_{on,i}, and considered the top 3,000 genes from this list. As a control, we chose the bottom 3,000 genes from the sorted list (Figure 4A). We found that genes with characteristics of rapid transcriptional bursting were associated with 'cell adhesion' (Figure 4B), consistent with previous reports that many tissuespecific cell adhesion molecules are expressed in mouse ES cells [37] and show celltocell variation in expression in mouse ES cell colonies [38]. Interestingly, the gene ontology category 'neural differentiation' was also enriched (Figure 4B), providing some support for previous studies that suggested that neural fate is chosen in a stochastic way [14, 39]. Conversely, the least varying set was enriched with genes associated with the maintenance of basic cellular function (Figure 4B), suggesting that eukaryotic cells have evolved to reduce the transcriptional noise of housekeeping genes for the phenotypic stability of basic cellular functions [14].
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
The steadystate solution of the chemical master equations can be written as a beta convolution of Poisson random variables (equation (3)). The mean and variance of the Poissonbeta distribution with these parameters are given by:
The squared coefficient of variation, η^{2}, and the Fano factor, φ, are given by:
A hierarchical Bayesian model
To estimate the parameters of the Poissonbeta distribution, we utilized a hierarchical Bayesian model, which can be described as:

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.
The full conditional distributions of s_{ i }, p_{ ij }, k_{on,i}, and k_{off,i}are nonstandard univariate. To sample the variables from their full conditional distributions, we use slice sampling [47]. For completeness, the full conditional distributions are given below:

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)$$
The log posterior probability, which can be used to monitor the convergence of the Gibbs sampler, is given by
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
To assess whether a set of observations generated from a Poissonbeta distribution follows a Poisson or negative binomial distribution, we used the parametric bootstrap for goodnessoffit testing [48]. We first generated n independent samples X_{1},..., X_{ n } from the Poissonbeta distribution with given parameters using the auxiliary variable representation. We then fitted these n simulated samples to the Poisson and negative binomial distributions using a maximum likelihood approach. The MATLAB function 'nbinfit' was used to compute the maximum likelihood estimates of the parameters of the negative binomial distribution. Based on the maximum likelihood estimates of the Poisson or negative binomial distribution, ${\theta}_{n}^{\mathsf{\text{dist}}}\left(\mathsf{\text{dist}}\in \left\{\mathsf{\text{Poisson}},\mathsf{\text{NB}}\right\}\right)$, we computed the KolmogorovSmirnov (KS) test statistic ${\mathsf{\text{KS}}}_{n}^{\mathsf{\text{dist}}}$ such that
where F_{ n } is the empirical distribution function for the n independent samples and ${F}_{{\theta}_{n}^{\mathsf{\text{dist}}}}$ is the cumulative distribution function of the Poisson or negative binomial distribution with the maximum likelihood estimates ${\theta}_{n}^{\mathsf{\text{dist}}}$. To evaluate the bootstrap P value, we repeated the following steps from k = 1 to k = B:

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).$$
Finally, the bootstrap P value is given by
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.
References
 1.
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.
 2.
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.
 3.
Wang Z, Gerstein M, Snyder M: RNASeq: a revolutionary tool for transcriptomics. Nature Reviews Genetics. 2009, 10: 5763. 10.1038/nrg2484.
 4.
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.
 5.
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.
 6.
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.
 7.
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.
 8.
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.
 9.
Tang F, Lao K, Surani MA: Development and applications of singlecell transcriptome analysis. Nature Methods. 2011, 8: S6S11. 10.1038/nchembio.740.
 10.
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.
 11.
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.
 12.
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.
 13.
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.
 14.
Eldar A, Elowitz MB: Functional roles for noise in genetic circuits. Nature. 2010, 467: 167173. 10.1038/nature09326.
 15.
Fuda NJ, Ardehali MB, Lis JT: Defining mechanisms that regulate RNA polymerase II transcription in vivo. Nature. 2010, 461: 186192.
 16.
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.
 17.
Peccoud J, Ycart B: Markovian modelling of gene product synthesis. Theoretical Population Biology. 1995, 48: 222234. 10.1006/tpbi.1995.1027.
 18.
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.
 19.
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.
 20.
Young RA: Control of the embryonic stem cell state. Cell. 2011, 144: 940954. 10.1016/j.cell.2011.01.032.
 21.
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.
 22.
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.
 23.
Silva J, Smith A: Capturing Pluripotency. Cell. 2008, 132: 532536. 10.1016/j.cell.2008.02.006.
 24.
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.
 25.
Munsky B, Neuert G, van Oudenaarden A: Using gene expression noise to understand gene regulation. Science. 2012, 336: 183187. 10.1126/science.1216379.
 26.
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.
 27.
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.
 28.
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.
 29.
Muller KE: Computing the confluent hypergeometric function, M(a, b, x). Numerishe Mathematik. 2001, 90: 179196. 10.1007/s002110100285.
 30.
Johnson NL, Kemp AW, Kotz S: Univariate discrete distributions. 2005, Wiley
 31.
Smiley MW, Proulx SR: Gene expression dynamics in randomly varying environments. Journal of Mathematical Biology. 2010, 61: 231251. 10.1007/s002850090298z.
 32.
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.
 33.
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.
 34.
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.
 35.
Suganuma T, Workman JL: Signals and combinatorial functions of histone modifications. Annual Review of Biochemistry. 2011, 80: 473499. 10.1146/annurevbiochem061809175347.
 36.
Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protocols. 2009, 4: 4457.
 37.
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.
 38.
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.
 39.
HemmatiBrivanlou A, Melton D: Vertebrate embryonic cells will become nerve cells unless told otherwise. Cell. 1997, 88: 1317. 10.1016/S00928674(00)81853X.
 40.
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.
 41.
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.
 42.
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.
 43.
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.
 44.
Miyanari Y, TorresPadilla M: Control of groundstate pluripotency by allelic regulation of Nanog. Nature. 2012, 483: 470473. 10.1038/nature10807.
 45.
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.
 46.
Anders S, Huber W: Differential expression analysis for sequence count data. Genome Biology. 2010, 11: R10610.1186/gb20101110r106.
 47.
Neal RM: Slice sampling. The Annals of Statistics. 2003, 31: 705767. 10.1214/aos/1056562461.
 48.
Best DJ, Rayner JCW: Goodness of fit for the Poisson distribution. Statistics & Probability Letters. 1999, 44: 259265. 10.1016/S01677152(99)000176.
 49.
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.
Acknowledgements
We thank Simon Anders, Wolfgang Huber, Anestis Touloumis, Nuno Fonseca and all members of the Marioni group for helpful comments.
Author information
Additional information
Competing interests
We declare no conflict of interests.
Authors' contributions
JKK performed all analysis and wrote all code. JKK and JCM wrote the manuscript and conceived the project. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Revised
Accepted
Published
DOI
Keywords
 gene regulation
 RNAseq
 singlecell
 statistics
 transcriptional burst