Inferring the kinetics of stochastic gene expression from single-cell RNA-sequencing data

Background

Genetically identical populations of cells grown in the same environmental condition show substantial variability in gene expression profiles. Although single-cell RNA-seq 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 single-cell RNA-seq data. By applying our model to a single-cell RNA-seq 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.

RNA-sequencing (RNA-seq) is a recently developed approach that allows an unbiased examination of the transcriptome to be performed using high-throughput DNA sequencing [1–3]. Compared to gene expression microarrays, the previous gold standard for genome-wide quantification of gene expression levels, RNA-seq 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 allele-specific expression to be interrogated [5, 6].

Until recently, most RNA-sequencing experiments began with a large population of cells (> 10⁵), 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, inter-cellular 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, single-cell 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 rate-limiting 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 rate-limiting 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 histone-modifying 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 pre-implantation 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 single-cell RNA-seq 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 single-cell RNA-sequencing 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.

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 rate-limiting 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 two-state 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:

Poisson-beta model. (A) Schematic of a two-state kinetic model for stochastic gene expression. (B) Heat map of the maximum P values of two goodness-of-fit tests for Poisson and negative binomial distributions. One thousand combinations of k _on and k _off were uniformly sampled from the log space by fixing s to 100. For each combination of the sampled parameters, 1,000 independent samples were generated from the Poisson-beta distribution to evaluate the fit of the data to the Poisson and negative binomial distributions using a bootstrap-based goodness-of-fit test. The colors represent minus log₁₀-transformed P values and the heat map is interpolated from the scattered data by using a Delaunay triangulation method. (C) Heat map of the Fano factor as a function of k _on and k _off with a fixed rate of transcription (s = 100). Along the black dashed line fixing the average number of mRNA molecules to 20, the four combinations of k _on and k _off give the varied level of the Fano factor and show different patterns of the variability of the number of mRNA molecules between cells. At point 1 with the highest Fano factor, the transitions between the two promoter states are slow, and the standardized expression level of a gene exhibits a U-shaped distribution, resulting in a bimodal distribution. At point 2, the transition to the inactive state is faster than the transition to the active state, and therefore the mRNA distribution has a long right tail resulting from occasional transcriptional bursts. As k _on and k _off increase at points 3 and 4, transitions between promoter states become fast, resulting in a Poisson-like distribution of the number of mRNA molecules with the Fano factor approaching 1. Note that this plot is similar to a recent figure generated by [25]. (D) Representative Poisson-beta distributions from four points in (C), which were computed with the auxiliary variable approach. (E) The corresponding beta distributions of p.

Full size image

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 Poisson-beta model

The parameters of the steady-state 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, ₁ F ₁, 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(x|s, k _on, k _off), which is known as the Poisson-beta distribution (PoBe) [30], takes the same form as the steady-state 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 two-state model, its stationary distribution takes the same form as that of p.

Given count measurements from RNA-sequencing 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 binomial-like 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 1B-E).

In practice, to ensure that the parameters are statistically identifiable, we suggest fitting three models (Poisson, negative binomial and Poisson-beta) to each gene before using a goodness-of-fit statistic to determine whether there is evidence that the parameters of the Poisson-beta 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 Poisson-beta model

Single-cell RNA-sequencing 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 Poisson-beta 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 single-cell RNA-seq 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 single-cell RNA-seq 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 population-based RNA-sequencing 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 spike-ins [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 Poisson-beta 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 non-identifiable parameters could be split into two broad categories (Additional file 1, Figure S3):

1. 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 Poisson-beta 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. 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 Poisson-beta model to data from only 12 independent ES cells.

To do this, we fitted the Poisson-beta 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 S4-S8), with a good agreement when 12 cells were considered (Additional file 1, Figure S6-S8). 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 S6-S8). 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 Poisson-beta 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 2H-I), 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).

Correlation of transcriptional kinetics with RNA polymerase II binding in mouse ES cells. In the left panel, burst size (or transcriptional efficiency) is plotted on the x-axis. In the middle panel, burst frequency is plotted on the x-axis. In the right panel, the average fraction of time that a gene spends in the inactive state is plotted on the x-axis. In all the panels, the following are plotted on the y-axis: the gene body activity (A)-(C), the gene promoter activity (D)-(F), and the pause index (G)-(I). Each point represents one identifiable gene with a normalized read count greater than 50 in at least one cell. ρ is the Spearman correlation coefficient.

Full size image

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 3A-I, 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 3J-L).

Correlation of transcriptional kinetics with histone modifications in mouse ES cells. (A)-(C) Box plots that compare burst size (A), burst frequency (B), and the average fraction of time that a gene spends in the inactive state (C) in the four groups. Given the annotated chromatin state of the two histone modifications by [33], we classified all expressed genes with a normalized read count greater than 50 in at least one cell that have chromatin state annotations into four groups: H3K4me3 only (n = 6,291), H3K27me3 only (n = 10), H3K4me3 + H3K27me3 (n = 630) and none (n = 492). The H3K4me3 group is significantly different from the others except the H3K27me3 group: P < 10^-16 for s _i /k _on,i , P < 10^-16 for k _on,i , P < 10^-16 for k _off,i /(k _off,i + k _on,i ), by the Mann-Whitney U-test. Due to the small number of samples of the H3K27me3 group, the H3K4me3 group is less significantly different from the H3K27me3 group: P = 0.0486 for s _i /k _on,i , P = 2.73 × 10^-5 for k _on,i , P = 8.12 × 10^-5 for k _off,i /(k _off,i + k _on,i ). In each box plot, the central red line indicates the median value, the top and bottom edges of the box are the 75th (q3) and 25th (q1) percentiles, and the ends of the whiskers denote q3 + 1.5(q3 - q1) and q1 - 1.5(q3 q1). (D)-(L) The profiles of H3K4me3, H3K27me3, and H3K36me3 ChIP-seq reads mapped near TSSs (Transcription Start Sites) are shown for genes with values of the three kinetic quantities in the top 10% (red), middle 10% (green), and bottom 10% (blue) of the relevant distribution. The y-axis is the total number of reads mapped to each position.

Full size image

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 tissue-specific cell adhesion molecules are expressed in mouse ES cells [37] and show cell-to-cell 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].

Functional enrichment of noisy genes showing distinctive bursting patterns of gene expression in mouse ES cells. (A) The mean of k _on,i and k _off,i (± S.E.) for the top and bottom 3,000 genes sorted by k _off,i / k _on,i . (S.E. = standard error). (B) The Benjamini corrected P values of the 18 GO terms enriched in the top or bottom 3,000 genes. GO, gene oncology.

Full size image

The Poisson-beta model provides a convenient statistical framework for modeling single-cell RNA-seq 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 Poisson-beta model are consistent with PolII binding and chromatin modifications using single-cell RNA-seq 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 single-cell RNA-seq studies have shown that such genes display less technical variability. However, as our understanding of the technical variability inherent to single-cell RNA-seq 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 rate-limiting step. While this assumption facilitates the derivation of a closed-form solution for the master equations, and thus the implementation of the Poisson-beta 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 4-thiouridine 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 RNA-seq [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 allele-specific expression when expression patterns are measured in single cells of a two-cell 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 allele-specific variation using the Poisson-beta model provides an interesting avenue for future research.

To summarize, as the single-cell 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 Poisson-beta model presented here, and other similar approaches, will be vital in maximizing the amount of biological insight that can be obtained from these data.

Properties of the steady-state solution of the master equation

The steady-state 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 Poisson-beta distribution with these parameters are given by:

The squared coefficient of variation, η², and the Fano factor, φ, are given by:

A hierarchical Bayesian model

To estimate the parameters of the Poisson-beta distribution, we utilized a hierarchical Bayesian model, which can be described as:

1. 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. 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. 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. 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. 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 top-level 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 non-standard 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. 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. 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. 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. 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).

Bootstrap-based goodness-of-fit test

To assess whether a set of observations generated from a Poisson-beta distribution follows a Poisson or negative binomial distribution, we used the parametric bootstrap for goodness-of-fit testing [48]. We first generated n independent samples X ₁,..., X _n from the Poisson-beta 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 Kolmogorov-Smirnov (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. 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. 2.

   Compute the maximum likelihood estimates ${\theta }_{n,k}^{*,\mathsf{\text{dist}}}$ from the bootstrap samples.

3. 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. 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 Poisson-beta 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 burn-in period.

Gene activity and pause index for GRO-Seq

To quantify PolII activity at the promoters and gene body regions, we defined three measures based on the number of mapped reads of GRO-Seq data [49]. First, gene body activity is defined as N/L where N is the number of GRO-Seq 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 Benjamini-corrected 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.

bp:

base pair

ES cell:

embryonic stem cell

GO:

gene ontology

ICM:

inner cell mass

PCR:

polymerase chain reaction

PoBe:

Poisson-beta

RNA FISH:

RNA fluorescence in situ hybridization

RNA PolII:

RNA polymerase II

RT:

reverse transcriptase.

We thank Simon Anders, Wolfgang Huber, Anestis Touloumis, Nuno Fonseca and all members of the Marioni group for helpful comments.

Authors and Affiliations

1. European Bioinformatics Institute (EMBL-EBI), Wellcome Trust Genome Sciences Campus, Hinxton, Cambridgeshire, CB10 1SD, UK

   Jong Kyoung Kim & John C Marioni

Corresponding author

Correspondence to John C Marioni.

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.

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.

Reprints and permissions