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

Abstract

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.

Figure 1

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.

Figure 2

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.

Figure 3

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.

Figure 4

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.