Single-cell RNA-seq reveals dynamic paracrine control of cellular variation (original) (raw)

Accession codes

Primary accessions

Gene Expression Omnibus

Data deposits

Data are deposited in GEO under accession number GSE48968.

References

1. Tay, S. et al. Single-cell NF-κB dynamics reveal digital activation and analogue information processing. Nature 466, 267–271 (2010)
   Article ADS CAS Google Scholar
2. Raj, A. & Van Oudenaarden, A. Single-molecule approaches to stochastic gene expression. Ann. Rev. Biophys. 38, 255–270 (2009)
   Article CAS Google Scholar
3. Slack, M. D., Martinez, E. D., Wu, L. F. & Altschuler, S. J. Characterizing heterogeneous cellular responses to perturbations. Proc. Natl Acad. Sci. USA 105, 19306–19311 (2008)
   Article ADS CAS Google Scholar
4. Sharma, S. V. et al. A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell 141, 69–80 (2010)
   Article CAS Google Scholar
5. Spencer, S. L., Gaudet, S., Albeck, J. G., Burke, J. M. & Sorger, P. K. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature 459, 428–432 (2009)
   Article ADS CAS Google Scholar
6. Taniguchi, Y. et al. Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science 329, 533–538 (2010)
   Article ADS CAS Google Scholar
7. Loewer, A. & Lahav, G. We are all individuals: causes and consequences of non-genetic heterogeneity in mammalian cells. Curr. Opin. Genet. Dev. 21, 753–758 (2011)
   Article CAS Google Scholar
8. Feinerman, O. et al. Single-cell quantification of IL-2 response by effector and regulatory T cells reveals critical plasticity in immune response. Mol. Syst. Biol. 6, 437–453 (2010)
   Article CAS Google Scholar
9. Veening, J.-W., Smits, W. K. & Kuipers, O. P. Bistability, epigenetics, and bet-hedging in bacteria. Ann. Rev. Microbiol. 62, 193–210 (2008)
   Article CAS Google Scholar
10. Fang, M., Xie, H., Dougan, S. K., Ploegh, H. & Van Oudenaarden, A. Stochastic cytokine expression induces mixed T helper cell states. PLoS Biol. 11, e1001618 (2013)
    Article CAS Google Scholar
11. Chalancon, G. et al. Interplay between gene expression noise and regulatory network architecture. Trends Genet. 28, 221–232 (2012)
    Article CAS Google Scholar
12. Bendall, S. C. et al. Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum. Science 332, 687–696 (2011)
    Article ADS CAS Google Scholar
13. Hashimshony, T., Wagner, F., Sher, N. & Yanai, I. CEL-Seq: single-cell RNA-Seq by multiplexed linear amplification. Cell Rep. 2, 666–673 (2012)
    Article CAS Google Scholar
14. Islam, S. et al. Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. Genome Res. 21, 1160–1167 (2011)
    Article CAS Google Scholar
15. Ramsköld, D. et al. Full-length mRNA-Seq from single-cell levels of RNA and individual circulating tumor cells. Nature Biotechnol. 30, 777–782 (2012)
    Article Google Scholar
16. Shalek, A. K. et al. Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature 498, 236–240 (2013)
    Article ADS CAS Google Scholar
17. Tang, F. et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nature Methods 6, 377–382 (2009)
    Article CAS Google Scholar
18. Amit, I. et al. Unbiased reconstruction of a mammalian transcriptional network mediating pathogen responses. Science 326, 257–263 (2009)
    Article ADS CAS Google Scholar
19. Wu, A. R. et al. Quantitative assessment of single-cell RNA-sequencing methods. Nature Methods 11, 41–46 (2014)
    Article CAS Google Scholar
20. Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. Proc. Natl Acad. Sci. USA 100, 9440–9445 (2003)
    Article ADS MathSciNet CAS Google Scholar
21. McDavid, A. et al. Data exploration, quality control and testing in single-cell qPCR-based gene expression experiments. Bioinformatics 29, 461–467 (2013)
    Article CAS Google Scholar
22. Islam, S. et al. Highly multiplexed and strand-specific single-cell RNA 5′ end sequencing. Nature Protocols 7, 823–828 (2012)
    Article Google Scholar
23. Marinov, G. K. et al. From single-cell to cell-pool transcriptomes: stochasticity in gene expression and RNA splicing. Genome Res. 24, 496–510 (2014)
    Article CAS Google Scholar
24. Garber, M. et al. A high-throughput chromatin immunoprecipitation approach reveals principles of dynamic gene regulation in mammals. Mol. Cell 47, 810–822 (2012)
    Article CAS Google Scholar
25. Weibrecht, I. et al. Visualising individual sequence-specific protein-DNA interactions in situ. New Biotechnol. 29, 589–598 (2012)
    Article CAS Google Scholar
26. Lee, T. K. et al. A noisy paracrine signal determines the cellular NFκB response to lipopolysaccharide. Sci. Signal. 2, ra65 (2009)
    PubMed PubMed Central Google Scholar
27. Rand, U. et al. Multi-layered stochasticity and paracrine signal propagation shape the type-I interferon response. Mol. Syst. Biol. 8, 584 (2012)
    Article Google Scholar
28. Zhao, M., Zhang, J., Phatnani, H., Scheu, S. & Maniatis, T. Stochastic expression of the interferon-β gene. PLoS Biol. 10, e1001249 (2012)
    Article CAS Google Scholar
29. Snijder, B. et al. Population context determines cell-to-cell variability in endocytosis and virus infection. Nature 461, 520–523 (2009)
    Article ADS CAS Google Scholar
30. Benveniste, E. N., Qin, H. & Type, I. Interferons as anti-inflammatory mediators. Sci. STKE 416, pe70 (2007)
    Google Scholar
31. Freedberg, I. M., Tomic-Canic, M., Komine, M. & Blumenberg, M. Keratins and the keratinocyte activation cycle. J. Invest. Dermatol. 116, 633–640 (2001)
    Article CAS Google Scholar
32. Rabani, M. et al. Metabolic labeling of RNA uncovers principles of RNA production and degradation dynamics in mammalian cells. Nature Biotechnol. 29, 436–442 (2011)
    Article CAS Google Scholar
33. Waters, C. M. & Bassler, B. L. Quorum sensing: cell-to-cell communication in bacteria. Ann. Rev. Cell Dev. Biol. 21, 319–346 (2005)
    Article CAS Google Scholar
34. Banchereau, J., Pascual, V. & Type, I. Interferon in systemic lupus erythematosus and other autoimmune diseases. Immunity 25, 383–392 (2006)
    Article CAS Google Scholar
35. Hall, J. C. & Rosen, A. Type I interferons: crucial participants in disease amplification in autoimmunity. Nature Rev. Rheumatol. 6, 40–49 (2010)
    Article CAS Google Scholar
36. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011)
    Article CAS Google Scholar

Download references

Acknowledgements

We thank B. Tilton, T. Rogers and M. Tam for assistance with cell sorting; E. Shefler, C. Guiducci, D. Thompson, and O. Rozenblatt-Rosen for project management and discussions and the Broad Genomics Platform for sequencing. We thank J. West, R. Lebofsky, A. Leyrat, M. Thu, M. Wong, W. Yorza, D. Toppani, M. Norris and B. Clerkson for contributions to C1 system development; B. Alvarado, M. Ray and L. Knuttson for assistance with C1 experiments; and M. Unger for discussions. Work was supported by an NIH Postdoctoral Fellowship (1F32HD075541-01, R.S.), an NIH grant (U54 AI057159, N.H.), an NIH New Innovator Award (DP2 OD002230, N.H.), an NIH CEGS (1P50HG006193-01, H.P., N.H., A.R.), NIH Pioneer Awards (5DP1OD003893-03 to H.P., DP1OD003958-01 to A.R.), the Broad Institute (H.P. and A.R.), HHMI (A.R.), the Klarman Cell Observatory at the Broad Institute (A.R.), an ISF-Broad Grant (N.F.), and the ERC (N.F.).

Author information

Author notes

1. Alex K. Shalek, Rahul Satija and Joe Shuga: These authors contributed equally to this work.

Authors and Affiliations

1. Department of Chemistry & Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, USA,
   Alex K. Shalek, Rona S. Gertner, Jellert T. Gaublomme, Ruihua Ding & Hongkun Park
2. Department of Physics, Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138, USA,
   Alex K. Shalek, Rona S. Gertner, Jellert T. Gaublomme, Ruihua Ding & Hongkun Park
3. Broad Institute of MIT and Harvard, 7 Cambridge Center, Cambridge, 02142, Massachusetts, USA
   Alex K. Shalek, Rahul Satija, John J. Trombetta, Dave Gennert, Diana Lu, Nir Yosef, Schraga Schwartz, Raktima Raychowdhury, Nir Hacohen, Hongkun Park & Aviv Regev
4. Fluidigm Corporation, 7000 Shoreline Court, Suite 100, South San Francisco, California 94080, USA,
   Joe Shuga, Peilin Chen, Brian Fowler, Suzanne Weaver, Jing Wang, Xiaohui Wang & Andrew P. May
5. School of Computer Science and Engineering, Hebrew University, 91904 Jerusalem, Israel,
   Nir Friedman
6. Center for Immunology and Inflammatory Diseases & Department of Medicine, Massachusetts General Hospital, Charlestown, 02129, Massachusetts, USA
   Nir Hacohen
7. Department of Biology, Howard Hughes Medical Institute, Massachusetts Institute of Technology, Cambridge, 02140, Massachusetts, USA
   Aviv Regev

Authors

1. Alex K. Shalek
   You can also search for this author inPubMed Google Scholar
2. Rahul Satija
   You can also search for this author inPubMed Google Scholar
3. Joe Shuga
   You can also search for this author inPubMed Google Scholar
4. John J. Trombetta
   You can also search for this author inPubMed Google Scholar
5. Dave Gennert
   You can also search for this author inPubMed Google Scholar
6. Diana Lu
   You can also search for this author inPubMed Google Scholar
7. Peilin Chen
   You can also search for this author inPubMed Google Scholar
8. Rona S. Gertner
   You can also search for this author inPubMed Google Scholar
9. Jellert T. Gaublomme
   You can also search for this author inPubMed Google Scholar
10. Nir Yosef
    You can also search for this author inPubMed Google Scholar
11. Schraga Schwartz
    You can also search for this author inPubMed Google Scholar
12. Brian Fowler
    You can also search for this author inPubMed Google Scholar
13. Suzanne Weaver
    You can also search for this author inPubMed Google Scholar
14. Jing Wang
    You can also search for this author inPubMed Google Scholar
15. Xiaohui Wang
    You can also search for this author inPubMed Google Scholar
16. Ruihua Ding
    You can also search for this author inPubMed Google Scholar
17. Raktima Raychowdhury
    You can also search for this author inPubMed Google Scholar
18. Nir Friedman
    You can also search for this author inPubMed Google Scholar
19. Nir Hacohen
    You can also search for this author inPubMed Google Scholar
20. Hongkun Park
    You can also search for this author inPubMed Google Scholar
21. Andrew P. May
    You can also search for this author inPubMed Google Scholar
22. Aviv Regev
    You can also search for this author inPubMed Google Scholar

Contributions

A.R., A.P.M., H.P., A.K.S., R.S. and J.S. conceived and designed the study. A.K.S., J.S., J.J.T., D.G., D.L., P.C., R.S.G., J.T.G., B.F., S.W., J.W., X.W., R.D. and R.R. performed experiments. R.S., A.K.S., S.S. and N.Y. performed computational analyses. R.S., A.K.S., J.S., N.F., H.P., A.P.M. and A.R. wrote the manuscript, with extensive input from all authors.

Corresponding authors

Correspondence toRahul Satija, Hongkun Park or Andrew P. May.

Ethics declarations

Competing interests

J.S., P.C., B.F., S.W., J.W., X.W. and A.P.M. declare competing financial interests as employees and/or stockholders in Fluidigm Corp.

Extended data figures and tables

Extended Data Figure 1 Single-cell RNA-seq of hundreds of dendritic cells.

a, Overview of experimental workflow. b, Shown are read densities for seven representative genes (two housekeeping genes (Rpl3 and Actb) and five immune response genes (Ifnb , Ifit1, Ccr7, Tnf, Marco)) across 60 single cells (blue) and one population control of 10,000 cells (grey; bulk population) after a 4 h LPS treatment. c, Distribution of failure scores for all single cells. Single cells with failure scores above 0.4 were discarded (see Supplementary Information). d–g, Comparisons of expression estimates for the average single cell and the bulk population. d, Scatter plots showing for each gene the relation between the average single-cell expression (y axis) and bulk population level expression (x axis) for each of four time points following LPS stimulation (1, 2, 4 and 6 h, left to right). e, The Pearson correlation coefficients (y axis) of each comparison, as in d, for each of the time points and stimuli presented in Fig. 1, as a function of the number of cells captured in the respective experiment (x axis). f, Scatter plots showing the residual (population-level expression minus the single cell average) in a LPS 1 h experiment (x axis) versus the residual in each of 3 other experiments (y axis, left to right): LPS 6 h, PIC 4 h and PAM 2 h. g, The Pearson correlation coefficient (y axis) between the bulk population level expression and the single-cell expression average when a different number of sub-sampled cells (x axis) is included in the single-cell average. h, i, Effects of Hoechst dye and periodic mixing on mRNA expression. h, Comparable expression levels after 4 h LPS with the addition of small amounts of Hoechst to aid in cell counting (x axis) and when no dye is used (y axis), when looking at all genes (left) or only immune response genes (right). i, Comparable expression levels after 4 h LPS with hourly mixing (x axis) or with no mixing (y axis), when looking at all genes (left) or only immune response elements (right). j, Core antiviral, peaked inflammatory, and sustained inflammatory module activation scores for a 0.1× LPS stimulation. Shown are violin plots of the scores (y axis) for the core antiviral (Supplementary Information, top), peaked inflammatory (Supplementary Information, middle), and sustained inflammatory modules (Supplementary Information, bottom) for each cell in (left to right): LPS 0 h, 1× (100 ng ml−1) LPS 4 h, and 0.1× (10 ng ml−1) LPS 4 h. k–n, PCA of stimulated dendritic cells. k, First two principal components (or PCs) from a PCA performed on the LPS stimulation time course. From top to bottom: unstimulated/LPS 0 h, LPS 1 h, LPS 2 h, LPS 4 h, LPS 6 h. l–n, PCAs (left) and the distributions of scores (right) for each of the first three PCs for samples collected after stimulation with LPS (l), PAM (m), or PIC (n), for 1 (yellow), 2 (blue), 4 (grey) and 6 (red) hours. A single PCA was performed for all cells in all three time courses.

Extended Data Figure 2 Effects of shallow read depth on expression estimates.

a, b, A million reads per cell are sufficient to estimate expression levels. a, Scatter plot for a single cell (from Shalek et al.16) showing the relation between expression estimates calculated using 30 M reads (x axis) or a sub-sample of 1 M reads (y axis). b, Scatter plots for six different dendritic cells stimulated for 4 h with LPS. Each plot shows the relation between expression estimates calculated using all reads (x axis; number of reads marked on axis label) or a sub-sample of 1 M reads (y axis). In all cases, R > 0.99. Note that although, in principle, no gene should be estimated as present only in a subsample but not the full data set, this does occur for a very small number of genes (for example, four genes in cell 3), representing a nuanced technical error in RNA-seq estimation. Consider two expressed genes, A and B, from distinct loci, but with a short stretch of sequence identity. At low sequencing depth, if reads only map to the shared region, estimation tools, such as RSEM36 (used here), can guess erroneously which gene is expressed, such that additional sequencing depth can ‘flip’ the assignment of an uncertain read from gene A to gene B. These cases are extremely rare, and have a negligible effect. c–e, A million reads per cell are sufficient to estimate μ, σ2, and α. Scatter plots showing the relation between α (c), μ (d), and σ2 (e) values estimated using 10 M reads per cell (on average; x axis) or a sub-sample of 1 M reads per cell (y axis) from RNA-seq libraries prepared from individual bone-marrow-derived dendritic cells stimulated for 4 h with LPS. c, For almost all genes, 1 M reads are sufficient to estimate α. For a very small fraction (<0.1%) of weakly expressed genes, estimates of α are improved with increased sequencing depth. For μ (d) and σ2 (e), estimates are plotted for all genes (left), only genes detected in more than 10 cells (middle), or only those genes detected in more than 30 cells (right).

Extended Data Figure 3 Technical and biological reproducibility.

a–d, Scatter plots showing the relationship between the average single-cell expression estimates in either of two technical replicates (LPS 2 h (a), LPS 4 h (b)) or two biological replicates (unstimulated/LPS 0 h (c), LPS 4 h (d)) for all genes (top), immune response genes (middle), or non-immune response genes (bottom). e, f, QQ plots (top) and MA plots (bottom) showing the similarity in expression estimates for the two LPS 2 h technical replicates (e) or the two LPS 4 h technical replicates (f). Plots are provided across all genes (left), non-immune response genes (middle), or immune response genes (right). g, h, QQ plots (top) and MA plots (bottom) showing the similarity in expression estimates for all cells (including cluster-disrupted cells) in the two LPS 0 h/unstimulated biological replicates (g) or the two LPS 4 h biological replicates (h) across either all genes (left), non-immune response genes (middle), or immune response genes (right). Note, that slight variations in the fraction of cluster-disrupted cells and activation state of one of the two 0 h samples results in mild deviations between immune response gene estimates in those biological replicates. i, j, PCA for the two LPS 4 h technical replicates. i, The first two principal components (PC1 and PC2, x and y axis, respectively) of a PCA from the two LPS 4 h stimulation technical replicates (blue: replicate 1; red: replicate 2). j, The distributions of scores for cells from each of the two technical replicates on each of the first five PCs (left to right: PC1 to PC5).

Extended Data Figure 4 Cluster disruption.

a, Single-cell expression distributions for SerpinB6b (a positive marker of cluster disruption) and Lyz1 (a negative marker of cluster disruption) at each time point (marked on top) after stimulation with LPS (all cells included, see Supplementary Information). Distributions are scaled to have the same maximum height. b, Difference in mRNA expression as measured by qRT–PCR (with a Gapdh control) between Lyz1 or SerpinB6b in cells pre-sorted before stimulation on the presence or absence of CD83 expression (CD83+ and CD83−, respectively), a known cell surface marker of cluster-disrupted cells (see Supplementary Information). Pre-sorted cells were then either unstimulated (blue) or stimulated (red) with LPS for 4 h. c, Expression of cluster-disruption markers does not change with stimulation. qRT–PCR showing the difference between Gapdh (control) and Lyz1 or SerpinB6b expression in cells pre-sorted on Cd83 either in the presence or absence of simulation with LPS. d, PCA showing the separation between maturing’ (blue dots) and cluster-disrupted (red dots) cells. e, Expression of cluster disruption markers for cells stimulated with LPS on-chip. For each cell (black dot) stimulated with LPS on-chip, shown are the expression levels (x axis) of SerpinB6b (cluster disruption cell marker, left) and Lyz1 (normal maturing cell marker, right) versus its antiviral score (y axis). With one exception, the cells are clearly maturing and not cluster-disrupted. Red shading: range of expression in cluster-disrupted cells.

Extended Data Figure 5 Time-dependent behaviours of single cells from different modules and stimuli.

a–c, For each of the three key modules: core antiviral, Id (a), peaked inflammatory, IIIc (b), and sustained inflammatory, IIId (c) shown are wave plots of all of its constituent genes in bone-marrow-derived dendritic cells stimulated with PAM (top), LPS (middle), or PIC (bottom) for 0, 1, 2, 4 and 6 h (left to right). x axis: expression level, ln(TPM + 1); y axis: genes; z axis: single-cell expression density (proportion of cells expressing at that level). Genes are ordered from lowest to highest average expression value at the 4 h LPS time point. d, Contributions of each module to measured variation. Significance of the contribution of modules Ia–Id and IIIa–IIId from Fig. 1 to the variation measured throughout the stimulation time course. Shown is the P value (Mann–Whitney test) of the tested association between each gene module and the first three PCs, calculated using a statistical resampling method (see Supplementary Information). Only the core antiviral, maturity, and peaked/sustained inflammatory clusters show statistically significant enrichments with the three PCs. e, Gene modules show coherent shifts in single-cell expression. Shown are heat maps of scaled α (left), μ (middle), and σ2 (right) values (colour bar, bottom) in each time course (LPS, PAM, PIC) for the genes in each of the three key modules (rows, modules marked on left). Heat maps are row-normalized across all three stimuli, with separate scalings for each of the three parameters, to highlight temporal dynamics. Genes are clustered as in Fig. 1. f, Dynamic changes in variation during stimulation for each module. For each module (rows) and stimulus (columns), shown are bar plots of the fraction of genes (y axis) with a significant change only in α (by a likelihood ratio test, P < 0.01, blue), only in μ (Wilcoxon test, P < 0.01, green), or in both (each test independently, light blue), at each transition (x axis), in different conditions (marked on top), separated by whether they increase or decrease during that transition. In each module and condition, the proportion is calculated out of the genes in the module that are significantly bimodal (by a likelihood ratio test) in at least one time point during the LPS response and are expressed in at least 10 cells in both conditions. Their number is marked on top of each bar; conditions with 3 or fewer genes changing are semi-transparent.

Extended Data Figure 6 Comparison of single-cell RNA-seq to RNA-FISH.

a–f, Single-cell mRNA expression distributions by RNA-FISH and single-cell RNA-seq. a, Representative images of genes analysed by RNA-FISH at 1 h and 4 h after LPS stimulation. b–f, mRNA expression distributions for the housekeeping gene B2m (b), the peaked inflammatory gene Cxcl1 (c), the core antiviral gene Ifit1 (d), the sustained inflammatory gene Il6 (e), and the peaked inflammatory gene Tnf (f) measured using either single-cell RNA-seq (top, black curves) or RNA-FISH (black histograms; no smoothing) in either unstimulated cells (LPS 0 h) or cells stimulated with LPS for 1, 2, 4 or 6 h. g–j, Determining the detection limit of single-cell RNA-seq by comparison to RNA-FISH. For each of 25 genes, we compared single-cell RNA-seq data (y axis, this study) to RNA-FISH data (x axis, from Shalek et al.) based on either their μ values (g, h) or α values (i, j), in which for RNA-FISH, expressing cells were defined for μ or α calculation at different thresholds (from left to right: at least 1, 4, 5, 10 or 20 copies per cell). g, i, Data from all 25 genes. h, j, Data after excluding probes from 5 cluster-disrupted gene markers (blue; Ccl22, Ccr7, Irf8, SerpinB6b, Stat4), which are less comparable since there are more cluster-disrupted cells in RNA-FISH experiments, and 2 low quality probes (grey; Pkm2, Fus) that showed very low expression in RNA-FISH, but had high expression in both single cell and population-level RNA-seq experiments. SPE (square-root of percent explained, top) represents the square-root of the total variance in the RNA-seq parameter explained by the RNA-FISH parameter, under the y = x model compared to a constant fit (that is, y = _y_–).

Extended Data Figure 7 Fitting gene expression distributions.

a, Flow chart of model fitting. Shown are the key steps in fitting the 3-parameter model. b, Examples of cases where fitting a multimodal distribution is required. Single-cell expression distributions for (top to bottom) Car13, Rgs1, Ms4a6c and Klf6 at (left to right) 1, 2, 4, and 6 h (marked on top) after stimulation with LPS. Distributions are scaled to have the same maximum height. Data: black lines; Bimodal fits: grey lines; Multimodal fits: blue lines. P values (colour-coded) calculated using a goodness-of-fit test (a low P value rejects the fit; see Supplementary Information). c–e, Reproducibility of gene-specific fitting of the undetected mode, when fitting a mix of two normal distributions to all data points, including those with ln(TPM + 1) < 1. c, d, Scatter plot showing the correlation between μ1 and μ2 estimates for the two LPS 4 h technical replicates (Supplementary Information), where μ1 and μ2 are the two component means (in decreasing order of magnitude) of the two mixed normal distributions. Estimates for μ2 correlate poorly between technical replicates, particularly when focusing on genes for which μ2 is greater than 1 (e), suggesting that the current data set does not support the use of this additional fit parameter. f, Robustness of α estimates to small deviations in the threshold. Scatter plots showing the correlation between α estimates determined when using a cut-off of ln(TPM + 1) = 1 (x axis) versus when using a cut-off of ln(TPM + 1) = 0.25 (y axis, left); 0.5 (y axis, middle) or 2 (y axis, right) for the LPS time course (top to bottom: 1 h, 2 h, 4 h and 6 h). g, Saturation curves for estimates of μ, σ2, and α. Box plots depicting the Pearson correlation coefficient between α (top), μ (middle), or σ2 (bottom) in two LPS 4 h technical replicates, as a function of the number of cells randomly drawn from each replicate (full details in Supplementary Information). Plots are shown for all genes (left), as well as those detected in more than 10 (middle) or 30 cells, (right) in both replicates (full data sets). h, i, Correcting for the relationship between mean expression and average detection. h, The probability of detecting a transcript (y axis) in a cell as a function of μ (x axis). Black, grey curves are two illustrative cells from the LPS 4 h time point. i, Differences in αMLE, a stringently-corrected MLE estimate of α (Supplementary Information), across the LPS time course. Shown are the box plots of αMLE values (y axis) for bimodally expressed genes (determined by a likelihood ratio test, Supplementary Information) at each time point (1, 2, 4, and 6 h) following LPS stimulation (x axis), as well as for the on-chip 4 h LPS stimulation, for each of the core antiviral (left), peaked inflammatory (middle) and sustained inflammatory (right) modules. Stars represent intervals where there is a significant difference in a parameter between two consecutive time points, as determined by a Wilcoxon rank sum test (single star: P < 10−2; double star: P < 10−5). j–l, Estimating an upper bound on α using a likelihood test. For each of three transcripts (Ifit1 (j); Rsad2 (k); and Cxcl1 (l)), shown are their expression distributions (red, left) and the matching likelihood function for a stringent upper estimate of α (blue dots, right), when considering a null model where expression is distributed in a lognormal fashion and any deviations are due to technical detection limits (Supplementary Information). Red vertical line: α MLE; black vertical line: nominal α. Vertical green bars signify the nominal estimation of α, representing the fraction of cells with detected expression of a transcript.

Extended Data Figure 8 Reproducibility of estimated μ, σ2 and α parameters.

a–f, Reproducibility of estimated μ, σ2 and α parameters between technical replicates. Scatter plots showing the relation between the estimated α (a), μ (b), and σ2 (c) values for the two unstimulated/LPS 0 h technical replicates. For μ (b) and σ2 (c), estimates are plotted for all genes (farthest on the left), as well as (left to right) for genes only detected in more than 10, 20, 30, 40 or 50 cells, respectively. d–f, show the same plots for the two LPS 4 h technical replicates. g, h, Reproducibility of estimated μ, σ2 and α parameters between biological replicates. Scatter plots showing the relation between the α (g), μ (h), and σ2 (i) values estimates for the two LPS 2 h biological replicates. For μ (h) and σ2 (i), estimates are plotted for all genes (farthest on the left), as well as (left to right) for genes only detected in more than 10, 20, 30, 40 or 50 cells, respectively. j–l, show the same plots for the two LPS 4 h biological replicates. m, n, Relationship between per-gene errors for μ, σ2 and α and the number of cells in which the gene’s expression is detected, or its bulk expression level. Scatter plots displaying the differences in the σ2 (left), μ (middle) and α (right) estimates for each gene between technical replicates for LPS 2 h (m) or LPS 4 h (n) (y axis) as a function of either the number of cells in which the transcript is detected (x axis, for μ and σ2), or the transcript’s bulk expression level (TPM, x axis, for α). Notably, σ2 (left) estimates saturate (denoted by a magenta line and shaded box) after ∼30 detectable events, whereas μ estimates saturate after ∼10. Dashed orange line: 95% confidence interval. o, p, Changes in μ, σ2 and α between time points are significant as compared to null models from both technical and biological replicates. Shown are the cumulative distribution functions (CDF) for shifts in μ (left), σ2 (middle), and α (right) between 2 h and 4 h (red CDF) for the core antiviral (top), peaked inflammatory (middle), and sustained inflammatory (bottom) modules compared to a null model (black CDF) derived from either technical (o) or biological (p) replicates (Supplementary Information). In the vast majority of cases, the changes between time points are significant, as assessed by a Kolmogorov–Smirnof (KS) test (P value of test in the upper left corner of each plot).

Extended Data Figure 9 Ifnb1 expression, production, and precocious cells.

a, b, Ifnb1 mRNA expression and the effect of IFN-β on variation. a, Single-cell expression distributions for the Ifnb1 transcript at each time point (top) after stimulation with PAM (top, green), LPS (middle, black), or PIC (bottom, magenta). Distributions are produced with the density function in R with default parameters, and scaled to have the same maximum density. b, For each of three modules defined in Fig. 1 (core antiviral, top; peaked inflammatory, middle; sustained inflammatory, bottom), shown are bar plots of the fraction of genes (y axis) with a significant change only in α (by a likelihood ratio test, P < 0.01, blue), only in μ (Wilcoxon test, P < 0.01, green), or in both (each test independently, light blue) between the 2 h LPS stimulation and the 2 h IFN-β stimulation separated by whether they increase or decrease during that transition. In each module and condition, the proportion is calculated out of the genes in the module that are significantly bimodal (by a likelihood ratio test) in at least time point during the LPS response and are expressed in at least 10 cells in both conditions. Their number is marked on top of each bar. c, d, Ifnb1 mRNA expression patterns and effect of cycloheximide. c, From top to bottom, population average Ifnb1 mRNA expression (top), single-cell average Ifnb1 mRNA expression (second to top), α (second to bottom) and μ (bottom) estimates for Ifnb1 for each stimulation condition in Fig. 1. Grey star at 6 h for PIC denotes uncertainty due to the small number of cells captured at that time point. d, Shown are box plots of the core antiviral scores (population level, see Supplementary Information) after a 4 h LPS stimulation either where cycloheximide was added at the time of stimulation (right, blue), or during a standard 4 h LPS control (left, green). Core antiviral expression is dramatically decreased by the addition of cycloheximide, suggesting that newly produced protein is needed to initiate the antiviral response. e, Relationship between core antiviral gene expression and Ifnb1 mRNA expression during the LPS, PAM and PIC stimulation time courses and in follow-up experiments. Shown are the expression of core antiviral genes (heat maps: rows, gene; columns, cells) for cells stimulated for 0, 1, 2, 4 or 6 h (left to right) with LPS (top), PAM (middle), or PIC (bottom). Beneath each heat map, grey bars depict the core antiviral score (middle panel, see Supplementary Information) and blue dots show Ifnb1 mRNA expression for each cell in every heat map. f–k, Identifying the precocious cells. f, Core antiviral scores for cells stimulated with LPS, PIC, or PAM. Shown are violin plots of the core antiviral module scores (Supplementary Information, y axis) for each cell from time course experiments (from left: 0, 1, 2, 4 and 6 h) of dendritic cells stimulated with LPS (top), PIC (middle) or PAM (bottom). Two precocious cells (yellow stars, top panel) have unusually high antiviral scores at 1 h LPS (yellow stars, top); similarly precocious cells can be seen in PIC at 1 h and 2 h (orange stars, middle) or in PAM at 2 h (turquoise stars, bottom). g, Precocious cells in all three responses are typical maturing cells. PCA showing the separation between maturing (blue dots) and cluster-disrupted (red dots) cells (top), as well as only maturing (middle) or only cluster-disrupted (bottom) cells (all as also shown in Extended Data Fig. 4d). The precocious cells from each of the responses are marked as stars (colours as in (f)), and all fall well within the maturing cells. h, Precocious cells in all three responses express Lyz1 and do not express SerpinB6b. Shown are mRNA expression distributions for SerpinB6b (cluster disruption cell marker, left) and Lyz1 (normal maturing cell marker, right) in LPS 1 h, PIC 1 h and 2 h, and PAM 1 h (top to bottom). The typical range for expression in cluster-disrupted cells is shaded in red. The precocious cells from each of the responses are marked as stars (colours as in (f)), and all fall well within the maturing cells. i, Normal quantile plots of the expression of genes from the core (cluster Id, left) and secondary (cluster Ic, right) antiviral clusters at 1 h LPS. The two precocious cells (yellow stars) express unusually high levels of core antiviral genes (left) but not of secondary genes (right). j, k, The precocious cells are only distinguished by the expression of ∼100 core antiviral genes. Shown are the distributions of scores for each of the first six PCs (right) for samples collected after stimulation with LPS for 1 h with (j) or without (k) the inclusion of core antiviral genes. Precocious cells (vertical red bars), normally distinguished by the third and fourth principle components (j), become indistinguishable from all other cells if the ∼ 100 core antiviral genes are excluded (k) before performing the PCA.

Extended Data Figure 10 Characterizing the precocious cells.

a, Core antiviral, peaked inflammatory, and sustained inflammatory module scores during the LPS time course and follow-up experiments. Shown are violin plots of the scores (y axis) for the core antiviral (Supplementary Information, top), peaked inflammatory (Supplementary Information, middle), and sustained inflammatory (Supplementary Information, bottom) modules for cells in each of the experiments (from left to right): LPS 0 h, LPS 1 h, LPS 2 h, LPS 2 h technical replicate 1, LPS 2 h technical replicate 2, LPS 4 h, LPS 4 h technical replicate 1, LPS 4 h technical replicate 2, LPS 4 h biological replicate, LPS 6 h, IFN-β 2 h, on-chip unstimulated, on-chip LPS 4 h, LPS 4 h with GolgiPlug at 0 h, LPS 4 h with GolgiPlug at 1 h, LPS 4 h with GolgiPlug at 2 h, LPS 4 h with _Ifnar_−/− dendritic cells, and LPS 4 h with _Stat1_−/− dendritic cells. Yellow stars: the two precocious cells at 1 h LPS. b, Changes in expression and variation during stimulation in the on-chip 4 h LPS stimulation. For genes in the (from top to bottom) core antiviral, maturity, peaked inflammatory and sustained inflammatory modules, shown are bar plots of the fraction of genes (y axis) with a significant change only in α (by a likelihood ratio test, P < 0.01, blue), only in μ (Wilcoxon test, P < 0.01, green), or in both (each test independently, light blue) between the 4 h on-chip LPS stimulation and the 4 h in-tube LPS stimulation separated by whether they increase or decrease during that transition. In each module and condition, the proportion is calculated out of the genes in the module that are significantly bimodal (by a likelihood ratio test) in at least one time point during the LPS response and are expressed in at least 10 cells in both conditions. Their number is marked on top of each bar. c–f, Bar plots, as in b, for a 4 h wild-type LPS stimulation (in-tube) and either a 4 h in-tube LPS stimulation where GolgiPlug was added 2 h after LPS (c), a 4 h LPS stimulation of _Ifnar_−/− dendritic cells (d), a 4 h LPS stimulation of _Stat1_−/− dendritic cells (e), or a 4 h LPS Stimulation of _Tnfr_−/− dendritic cells (f). g, Genetic perturbations alter expression and variation in different inflammatory and antiviral modules. Shown is the expression of the genes (rows) in, from top to bottom: the core antiviral (Id), maturity (IIIb), peaked inflammatory (IIIc), and sustained inflammatory (IIId) modules in single cells (columns) in, from left to right: the in-tube, on-chip, _Ifnar1_−/−, _Stat1_−/−, and _Tnfr_−/− conditions. Yellow/purple colour scale: scaled expression values (_z_-scores). h, Scores of the peaked inflammatory module for _Ifnar_−/− dendritic cells with recombinant cytokines. Shown are box plots of the peaked inflammatory module scores (Supplementary Information, y axis) for three population-level replicates of a 4 h LPS stimulation of _Ifnar_−/− dendritic cells to which a recombinant cytokine (x axis) has been added at 2 h after stimulation. Notably, adding IL-10 significantly reduces the peaked inflammatory module. i, Stat1 knockout affects expression and variation of peaked inflammatory genes. Shown are expression distributions for five peaked inflammatory genes after 4 h of LPS stimulation in each of three conditions: in-tube stimulation of wild-type dendritic cells (control; left), on-chip stimulation of wild-type cells (no cell-to-cell signalling; middle), and a stimulation of dendritic cells from _Stat1_−/− mice (performed in-tube; right). j, k, Population-level paracrine signalling enhances and coordinates the core antiviral response while dampening and desynchronizing the peaked inflammatory ones. j, Gene network model showing how positive IFN-β signalling induces the antiviral response and reduces its heterogeneity, while simultaneously activating negative paracrine feedback, possibly including IL-10, which dampens the peaked inflammatory cluster and increases its heterogeneity. k, Cell population model showing how positive and negative paracrine signalling alter antiviral (magenta) and inflammatory (green) gene expression variability across cells. Grey denotes no expression.

Supplementary information

Supplementary Information

This file contains Supplementary Notes, Supplementary Methods and Supplementary References. (PDF 855 kb)

Supplementary Table 1

This table contains a summary of single cell RNA-Seq experiments. (XLS 34 kb)

Supplementary Table 2

This table contains sequencing metrics for single cell and population RNA-seq libraries. (XLS 150 kb)

Supplementary Table 3

This table shows cluster assignments and GO enrichments for the 813 genes that were induced by at least two-fold (compared to unstimulated cells) at a population level at any time point during the LPS time course. (XLS 671 kb)

Supplementary Table 4

This table contains ingenuity analysis of the genes differentially regulated between a wildtype 4 hour LPS stimulation and either the "on-chip" stimulation, an Ifnar knockout, or a Stat1 knockout. (XLS 1657 kb)

PowerPoint slides

Rights and permissions

About this article

Cite this article

Shalek, A., Satija, R., Shuga, J. et al. Single-cell RNA-seq reveals dynamic paracrine control of cellular variation.Nature 510, 363–369 (2014). https://doi.org/10.1038/nature13437

Download citation

• Received: 16 June 2013
• Accepted: 02 May 2014
• Published: 11 June 2014
• Issue Date: 19 June 2014
• DOI: https://doi.org/10.1038/nature13437