Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer

Abstract

Modulation of chromatin structure via histone modification is a major epigenetic mechanism and regulator of gene expression. However, the contribution of chromatin features to tumor heterogeneity and evolution remains unknown. Here we describe a high-throughput droplet microfluidics platform to profile chromatin landscapes of thousands of cells at single-cell resolution. Using patient-derived xenograft models of acquired resistance to chemotherapy and targeted therapy in breast cancer, we found that a subset of cells within untreated drug-sensitive tumors share a common chromatin signature with resistant cells, undetectable using bulk approaches. These cells, and cells from the resistant tumors, have lost chromatin marks—H3K27me3, which is associated with stable transcriptional repression—for genes known to promote resistance to treatment. This single-cell chromatin immunoprecipitation followed by sequencing approach paves the way to study the role of chromatin heterogeneity, not just in cancer but in other diseases and healthy systems, notably during cellular differentiation and development.

This is a preview of subscription content, access via your institution

Access options

Rent or buy this article

Prices vary by article type

from$1.95

to$39.95

Prices may be subject to local taxes which are calculated during checkout

Fig. 1: Reconstructing cell-type-specific chromatin states from single-cell ChIP-seq profiles.
Fig. 2: Single-cell ChIP-seq profiling of mouse stromal cells reveals cell identities from H3K27me3 chromatin landscapes.
Fig. 3: Sensitive and drug-resistant specific H3K27me3 chromatin landscapes in PDX model of triple-negative breast cancer treated with capecitabine.
Fig. 4: A fraction of cells from sensitive tumor share H3K27me3 chromatin features with resistant cells in a model of luminal estrogen receptor-positive PDX treated with tamoxifen.

Similar content being viewed by others

Data availability

All sequencing files and processed count matrices were deposited with the Gene Expression Omnibus under accession number GSE117309.

Code availability

Codes are available from the GitHub repository (https://github.com/vallotlab/scChIPseq).

References

  1. Roth, A. et al. Clonal genotype and population structure inference from single-cell tumor sequencing. Nat. Methods 13, 573–576 (2016).

    Article  CAS  Google Scholar 

  2. Nik-Zainal, S. et al. The life history of 21 breast cancers. Cell 149, 994–1007 (2012).

    Article  CAS  Google Scholar 

  3. McGranahan, N. & Swanton, C. Clonal heterogeneity and tumor evolution: past, present, and the future. Cell 168, 613–628 (2017).

    Article  CAS  Google Scholar 

  4. Dagogo-Jack, I. & Shaw, A. T. Tumour heterogeneity and resistance to cancer therapies. Nat. Rev. Clin. Oncol. 15, 81–94 (2018).

    Article  CAS  Google Scholar 

  5. Salgia, R. & Kulkarni, P. The genetic/non-genetic duality of drug ‘resistance’ in cancer. Trends Cancer 4, 110–118 (2018).

    Article  CAS  Google Scholar 

  6. Gardner, E. E. et al. Chemosensitive relapse in small cell lung cancer proceeds through an EZH2-SLFN11 axis. Cancer Cell 31, 286–299 (2017).

    Article  CAS  Google Scholar 

  7. Shaffer, S. M. et al. Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature 546, 431–435 (2017).

    Article  CAS  Google Scholar 

  8. Rathert, P. et al. Transcriptional plasticity promotes primary and acquired resistance to BET inhibition. Nature 525, 543–547 (2015).

    Article  CAS  Google Scholar 

  9. Kim, C. et al. Chemoresistance evolution in triple-negative breast cancer delineated by single-cell sequencing. Cell 173, 879–893.e13 (2018).

    Article  CAS  Google Scholar 

  10. Horning, A. M. et al. Single-cell RNA-seq reveals a subpopulation of prostate cancer cells with enhanced cell-cycle-related transcription and attenuated androgen response. Cancer Res. 78, 853–864 (2018).

    Article  CAS  Google Scholar 

  11. Mazor, T. et al. DNA methylation and somatic mutations converge on the cell cycle and define similar evolutionary histories in brain tumors. Cancer Cell 28, 307–317 (2015).

    Article  CAS  Google Scholar 

  12. Aryee, M. J. et al. DNA methylation alterations exhibit intraindividual stability and interindividual heterogeneity in prostate cancer metastases. Sci. Transl. Med. 5, 169ra10 (2013).

    Article  Google Scholar 

  13. Clark, S. J. et al. Genome-wide base-resolution mapping of DNA methylation in single cells using single-cell bisulfite sequencing (scBS-seq). Nat. Protoc. 12, 534–547 (2017).

    Article  CAS  Google Scholar 

  14. Angermueller, C. et al. Parallel single-cell sequencing links transcriptional and epigenetic heterogeneity. Nat. Methods 13, 229–232 (2016).

    Article  CAS  Google Scholar 

  15. Li, L. et al. Single-cell multi-omics sequencing of human early embryos. Nat. Cell Biol. 20, 847–858 (2018).

    Article  CAS  Google Scholar 

  16. Cusanovich, D. A. et al. The cis-regulatory dynamics of embryonic development at single-cell resolution. Nature 555, 538–542 (2018).

    Article  CAS  Google Scholar 

  17. Kundaje, A. et al. Integrative analysis of 111 reference human epigenomes. Nature 518, 317–330 (2015).

    Article  CAS  Google Scholar 

  18. Ernst, J. et al. Mapping and analysis of chromatin state dynamics in nine human cell types. Nature 473, 43–49 (2011).

    Article  CAS  Google Scholar 

  19. Rotem, A. et al. Single-cell ChIP-seq reveals cell subpopulations defined by chromatin state. Nat. Biotechnol. 33, 1165–1172 (2015).

    Article  CAS  Google Scholar 

  20. Marangoni, E. et al. Capecitabine efficacy is correlated with TYMP and RB1 expression in PDX established from triple-negative breast cancers. Clin. Cancer Res. 24, 2605–2615 (2018).

    Article  CAS  Google Scholar 

  21. Lederer, M., Bley, N., Schleifer, C. & Hüttelmaierr, S. The role of the oncofetal IGF2 mRNA-binding protein 3 (IGF2BP3) in cancer. Semin. Cancer Biol. 29, 3–12 (2014).

    Article  CAS  Google Scholar 

  22. Zheng, X. et al. Epithelial-to-mesenchymal transition is dispensable for metastasis but induces chemoresistance in pancreatic cancer. Nature 527, 525–530 (2015).

    Article  CAS  Google Scholar 

  23. Fischer, K. R. et al. Epithelial-to-mesenchymal transition is not required for lung metastasis but contributes to chemoresistance. Nature 527, 472–476 (2015).

    Article  CAS  Google Scholar 

  24. Cottu, P. et al. Acquired resistance to endocrine treatments is associated with tumor-specific molecular changes in patient-derived luminal breast cancer xenografts. Clin. Cancer Res. 20, 4314–4325 (2014).

    Article  CAS  Google Scholar 

  25. Massarweh, S. et al. Tamoxifen resistance in breast tumors is driven by growth factor receptor signaling with repression of classic estrogen receptor genomic function. Cancer Res. 68, 826–833 (2008).

    Article  CAS  Google Scholar 

  26. Ciupek, A. et al. Androgen receptor promotes tamoxifen agonist activity by activation of EGFR in ERα-positive breast cancer. Breast Cancer Res. Treat. 154, 225–237 (2015).

    Article  CAS  Google Scholar 

  27. Denduluri, S. K. et al. Insulin-like growth factor (IGF) signaling in tumorigenesis and the development of cancer drug resistance. Genes Dis. 2, 13–25 (2015).

    Article  CAS  Google Scholar 

  28. Nguyen, V. T. et al. Differential epigenetic reprogramming in response to specific endocrine therapies promotes cholesterol biosynthesis and cellular invasion. Nat. Commun. 6, 10044 (2015).

    Article  CAS  Google Scholar 

  29. Magnani, L. et al. Genome-wide reprogramming of the chromatin landscape underlies endocrine therapy resistance in breast cancer. Proc. Natl Acad. Sci. USA 110, E1490–E1499 (2013).

    Article  CAS  Google Scholar 

  30. Cottu, P. et al. Modeling of response to endocrine therapy in a panel of human luminal breast cancer xenografts. Breast Cancer Res. Treat. 133, 595–606 (2012).

    Article  CAS  Google Scholar 

  31. Marangoni, E. et al. A new model of patient tumor-derived breast cancer xenografts for preclinical assays. Clin. Cancer Res. 13, 3989–3998 (2007).

    Article  CAS  Google Scholar 

  32. Petit, V. et al. Optimization of tumor xenograft dissociation for the profiling of cell surface markers and nutrient transporters. Lab. Invest. 93, 611–621 (2013).

    Article  CAS  Google Scholar 

  33. Duffy, D. C., McDonald, J. C., Schueller, O. J. & Whitesides, G. M. Rapid prototyping of microfluidic systems in poly(dimethylsiloxane). Anal. Chem. 70, 4974–4984 (1998).

    Article  CAS  Google Scholar 

  34. Mazutis, L. et al. Single-cell analysis and sorting using droplet-based microfluidics. Nat. Protoc. 8, 870–891 (2013).

    Article  CAS  Google Scholar 

  35. Siegel, A. C., Bruzewicz, D. A., Weibel, D. B. & Whitesides, G. M. Microsolidics: fabrication of three-dimensional metallic microstructures in poly(dimethylsiloxane). Adv. Mater. 19, 727–733 (2007).

    Article  CAS  Google Scholar 

  36. Clausell-Tormos, J. et al. Droplet-based microfluidic platforms for the encapsulation and screening of mammalian cells and multicellular organisms. Chem. Biol. 15, 427–437 (2008).

    Article  CAS  Google Scholar 

  37. Anna, S. L., Bontoux, N. & Stone, H. A. Formation of dispersions using “flow focusing” in microchannels. Appl. Phys. Lett. 82, 364–366 (2003).

    Article  CAS  Google Scholar 

  38. Zilionis, R. et al. Single-cell barcoding and sequencing using droplet microfluidics. Nat. Protoc. 12, 44–73 (2017).

    Article  CAS  Google Scholar 

  39. Klein, A. M. et al. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell 161, 1187–1201 (2015).

    Article  CAS  Google Scholar 

  40. Abate, A. R., Chen, C. H., Agresti, J. J. & Weitz, D. A. Beating Poisson encapsulation statistics using close-packed ordering. Lab Chip 9, 2628–2631 (2009).

    Article  CAS  Google Scholar 

  41. Chabert, M., Dorfman, K. D. & Viovy, J. L. Droplet fusion by alternating current (AC) field electrocoalescence in microchannels. Electrophoresis 26, 3706–3715 (2005).

    Article  CAS  Google Scholar 

  42. Vallot, C. et al. Erosion of X chromosome inactivation in human pluripotent cells initiates with XACT coating and depends on a specific heterochromatin landscape. Cell Stem Cell 16, 533–546 (2015).

    Article  CAS  Google Scholar 

  43. Liu, Y., Popp, B. & Schmidt, B. CUSHAW3: sensitive and accurate base-space and color-space short-read alignment with hybrid seeding. PLoS ONE 9, e86869 (2014).

    Article  Google Scholar 

  44. Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25 (2009).

    Article  Google Scholar 

  45. Wilkerson, M. D. & Hayes, D. N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573 (2010).

    Article  CAS  Google Scholar 

  46. van der Maaten, L. & Hinton, G. Visualizing data using t-SNE. J. Mach. Learn. Res. 9, 2579–2605 (2008).

    Google Scholar 

  47. Phanstiel, D. H., Boyle, A. P., Araya, C. L. & Snyder, M. P. Sushi.R: flexible, quantitative and integrative genomic visualizations for publication-quality multi-panel figures. Bioinformatics 30, 2808–2810 (2014).

    Article  CAS  Google Scholar 

  48. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B 57, 289–300 (1995).

    Google Scholar 

  49. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).

    Article  CAS  Google Scholar 

  50. Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005).

    Article  CAS  Google Scholar 

  51. Khandelwal, G. et al. Next-generation sequencing analysis and algorithms for PDX and CDX models. Mol. Cancer Res. 15, 1012–1016 (2017).

    Article  CAS  Google Scholar 

  52. McCarthy, D. J., Campbell, K. R., Lun, A. T. & Wills, Q. F. Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R. Bioinformatics 33, 1179–1186 (2017).

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Barron, M. & Li, J. Identifying and removing the cell-cycle effect from single-cell RNA-sequencing data. Sci. Rep. 6, 33892 (2016).

    Article  CAS  Google Scholar 

  54. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).

    Article  CAS  Google Scholar 

  55. Lai, D., Ha, G. & Shah, S. HMMcopy: copy number prediction with correction for GC and mappability bias for HTS data. R package v.1.8.0 (2016); http://bioconductor.org/packages/release/bioc/html/HMMcopy.html

Download references

Acknowledgements

We thank members of the Griffiths and Vallot laboratories and HiFiBiO for helpful discussions. We thank B. E. Bernstein (Massachusetts General Hospital, Broad Institute) and R. Nicol (Whitehead Institute, Broad Institute) for advice on technology development, J. Waterfall (Institut Curie) and B.E.B. for critical discussion on the manuscript and the Institut Pierre-Gilles de Gennes for the use of the clean room facilities. We also thank patients for their approval to use their tumor samples. This work was supported by the Investissements d’Avenir program via Agence Nationale de la Recherche grant nos. ANR-10-NANO-02, ANR-10-IDEX-0001-02 PSL, ANR-10-LABX-31, ANR-10-EQPX-34, ANR-10-EQPX-03, ANR-14-CE10-0013_CellectChIP and ANR-10-INBS-09-08 (France Génomique Consortium). K.G. was supported by a Cifre PhD funded by ANRT and HiFiBiO, fellowship no. 2015/0159. C.V. was supported by the ATIP-Avenir program, by the Plan Cancer program and by the Institut Curie SiRIC program, SiRIC grant nos. INCa-DGOS- 4654 and INCa-DGOS-Inserm_12554. The scRNA-seq experiments were supported by the PSL University grant ‘Single-Cell Sequencing’. Next-generation sequencing was performed by the ICGex platform of the Institut Curie and the iGenSeq platform of the Institut du Cerveau et de la Moelle épinière.

Author information

Authors and Affiliations

Authors

Contributions

The scChIP-seq experimental method was conceived by K.G., C.B. and A.G., developed by K.G. and A.G. for the cellular and molecular biology and developed by K.G., A.D.G and A.G. for the microfluidics. The scChIP-seq analytical methods were conceived and developed by K.G. and C.V.; K.G., C.B., A.D.G., C.V. and A.G. designed the experiments. K.G., A.D., F.N., A.D., A.P. and O.F. conducted the experiments. The PDX experiments were designed and interpreted by J.M., F.R., A.D. and E.M.; S.L. performed the sequencing experiments. Y.P. performed the barcoded bead preparation and quality controls. M.R. contributed to the microfluidics chip design. K.G., A.W. and C.V. performed the data analysis. K.G., A.D.G., C.V. and A.G. wrote the manuscript with input from all authors.

Corresponding authors

Correspondence to Andrew D. Griffiths, Céline Vallot or Annabelle Gérard.

Ethics declarations

Competing interests

Patents have been filed on some aspects of this work. The inventors may receive payments related to the exploitation of these patents under their employer’s rewards to inventor scheme. A.D.G. and C.B. are cofounders of HiFiBiO.

Additional information

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Integrated supplementary information

Supplementary Fig. 1 Overview of the single-cell ChIP-seq procedure.

(a) The single-cell ChIP-seq procedure is divided into 3 workflows including a droplet-microfluidic workflow, a ChIP-seq workflow and an analytical workflow. Key steps of each workflow are indicated in the boxes. (b) Single-cell ChIP-seq timeline: the entire process from the cell harvesting to the NGS takes 5 days. IP: immunoprecipitation; QC: Quality Controls; NGS: Next Generation Sequencing.

Supplementary Fig. 2 Design of microfluidic devices.

(a) Device used to compartmentalize cells in 45 pl droplets. (b) Device used to produce 9 pl droplets to make hydrogel beads. (c) Device used for compartmentalization of hydrogel beads in 100 pl droplets. (d) Device used to merge nucleosome-containing droplets with hydrogel barcoded bead- containing droplets. Potential and ground electrodes are indicated with the ‘+’ and the ‘-’ marks, respectively. Scale bars correspond to 2 mm.

Supplementary Fig. 3 Barcoded bead production and quality control.

(a) Beads were produced in a microfluidic device with 2-inlets by dispersing a mixture comprising PolyEthylene Glycol Di-Acrydrite (PEG-DA), Streptavidin Acrylamide and the photo-initiator. Flow rates were adjusted to produce 9 pl droplets, and immediately exposed to UV light for polymerization of the hydrogel network (see Methods). Scale bar corresponds to 25 μm. (b) Split- and-pool synthesis principle for the addition of successive indexes. (c) Barcodes were synthesized by successive ligation of double-stranded indexes containing 5′ overhang of 4 base pairs by three rounds of split-and-pool synthesis using 96 Index 1, 96 Index 2 and 96 Index 3. Barcodes were flanked at one end by common sequences comprising a 1⁄2 Pac1 restriction site, a T7 promoter and the Illumina Read #2 sequencing primer, which were bound to the beads via a photocleavable linker (PC-linker). A 3′ C3-spacer was added to the 3′end of the photocleaved site for directed ligation to the other end of the barcode comprising a second common sequence with the 1⁄2 Pac1 restriction site ligated to the index 3. (d) Barcodes that failed in one of the three split-pool rounds were completed with a ‘block’ oligonucleotide comprising a 5′ C3-spacer and a 3′ Inverted ddT to prevent ligation. (e) Single-bead sequencing results showing the fraction of the first two most abundant barcodes of 16 beads. In average, 97.7% of the barcodes corresponded to the same sequence and the second most abundant barcode represented 0.17% of all sequences.

Supplementary Fig. 4 Sequencing library preparation.

(a) Enriched barcoded nucleosomes were linearly amplified by in vitro transcription. The amplified RNAs were reverse transcribed into cDNA by random priming, appending a reverse complement of Illumina Read #1 sequencing primer. The cDNAs were amplified by PCR, appending an Illumina P7 and P5 sequences. (b) Schematic of the final sequencing product with size in bp of each element constituting the sequence. (c) Electropherogram showing the size distribution of the final sequencing library post agarose gel purification. The smear ranges from 300 bp to 700 bp and corresponds to barcoded nucleosomes (profile obtained by Tapestation [Agilent]). (d) Single-cell ChIP-seq libraries were sequenced as follows: 50 bp were assigned to read the nucleosomal sequence and 100 bp were assigned to read the barcode.

Supplementary Fig. 5 Live monitoring of the fusion between droplets with cells and with barcoded hydrogel beads on the microfluidic system.

(a) Time traces of 5 droplets showing possible fusion events occurring in all scChIP-seq experiments: droplet #1 contained a cell but no bead; droplet #2 and droplet #5 contained one bead but no cell; droplet #3 is empty; droplet #4 contained one cell plus one bead, resulting in nucleosomes barcoding. Droplets were scanned after fusion as they crossed a laser beam and their fluorescence intensity analyzed in real time. Sulforhodamine B (orange fluorescent) was used as common marker (drop-code) and Dye-405 (violet fluorescent) was used as marker of the nucleosome-containing droplets (drop-code ‘Cell’). Cells and beads were labeled with Calcein AM cellpermeant dye (green fluorescent) and Biotin-Cy5 (red fluorescent, ‘Bead-code’), respectively. (b) Scatter plot showing cell fluorescence intensity (green, ‘Cell-code’) versus bead fluorescence intensity (red, ‘Bead-code’) in each droplet allowing precise counting of the number of cells co-encapsulated with a barcoded bead. Droplets from panel (a) are indicated as examples of the different merged droplet populations. (c) Barplot showing for two scChIP-seq experiments targeting distinct histone modifications (H3K27me3 and H3K4me3), the total number of encapsulated cells detected by fluorescence on the microfluidic station, the number of cells coencapsulated with a barcoded hydrogel bead detected by fluorescence on the microfluidic station and the number of cells (that is unique barcodes) identified from the sequencing data.

Supplementary Fig. 6 Human and mouse cells mixture confirms single-cell resolution.

(a) Scatter plot of number of reads per barcode aligning to the mouse versus human reference genome showing that 96.5% of the barcodes are specific to one species (at least 95% of the reads with the same barcode mapping to one of the two species). The percentage of mouse (26.4%), human (70.1%) and mixed (3.5%) species was close to expected values based on a Poisson distribution of cells in droplets with an average number of cells per droplets l of 0.1 (32.6%, 65.2% and 2.2%, respectively). (b) Barplot showing the number of barcodes (corresponding to the number of cells) identified for each species (blue, red and black bars) in comparison to the expected number of cells counted on the microfluidic station (gray bars; 3,000 in total from a mixture comprising 1/3 mouse cells and 2/3 of human cells).

Supplementary Fig. 7 H3K4me3 and H3K27me3 single-cell ChIP-seq datasets for in vitro experiments.

(a) Histograms of the distribution of scChIP-seq raw and unique sequencing reads per barcode (that is per cell) in H3K4me3 and H3K27me3 single-cell ChIP-seq datasets. (b) Density scatterplots representing log2 cumulative counts between three independent fractions of the same emulsion of B cells collected and processed in parallel to generate H3K4me3 single-cell ChIP-seq datasets. Correlation between replicates is calculated based on the cumulative count per million reads in 5 kb genomic bins across single-cells (n = 642,098 genomic bins). Pearson’s correlation scores and P values are computed genome-wide. (c) Density scatter plot for two biological replicates corresponding to two emulsions of B cells collected from different cell culture flasks and processed with different batches of barcoded hydrogel beads to generate H3K27me3 single-cell ChIP-seq datasets. Correlation between replicates is calculated based on the cumulative count per million reads in 50 kb genomic bins across single-cells (n = 64,455 genomic bins). Pearson’s correlation score and P value are computed genome-wide. (d) t-SNE plots representing H3K27me3 scChIP-seq data from two biological replicates (n = 1,813 and 1,451 single-cells in replicate #1 and #2 respectively), colored according to the batch of origin (left) or consensus clustering results (right), arguing against a batch effect for cell population clustering. (e) Left panel: hierarchical clustering and corresponding heatmap of cell-to-cell Pearson’s correlation scores for H3K27me3 scChIP-seq data from a mixed population of 1:1 human B cells and T cells together (n = 1,453 single-cells) with separately barcoded B cells (n = 922 single-cells) and T cells (n = 879 single-cells). The unique read count, the batch of origin and the consensus clustering result is indicated above the heatmap. Right panel: Corresponding t-SNE plot, points from the mixed population are colored in gray, while points from the separately barcoded B cells and T cells are colored according to the cell type-specific barcode sequence. (f) Venn diagrams comparing for T and B cells datasets the H3K4me3 peaks detected by single-cell and bulk approaches.

Supplementary Fig. 8 Single-cell ChIP-seq data distinguish human T cells (Jurkat) from human B cells (Ramos).

(a) Consensus clustering matrix for H3K4me3 (top panel) and H3K27me3 (bottom panel) scChIP-seq datasets. Consensus score ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). (b) Volcano plots representing adjusted P values (two-sided Wilcoxon signed-rank test) versus fold-changes for differential analysis comparing chromatin features between B cells and T cells (thresholds of 0.01 for adjusted P value and 1 for |log2FC|) for H3K4me3 (top panel, n = 1,772 enriched and 926 depleted regions) and H3K27me3 (bottom panel, n = 2,126 enriched and 7,327 depleted regions) scChIP-seq datasets. (c) Barplot displaying the -log10 of adjusted P values from pathway analysis in H3K4me3 scChIP-seq dataset. The top 10 significant gene sets are indicated below the barplot.

Supplementary Fig. 9 Residual triple-negative breast tumor model of acquired resistance to Capecitabine.

(a) n = 6 mice with HBCx-95 derivatives were treated with Capecitabine for 6 weeks. Graphs represent relative tumor volumes (RTV, mm3) with time. Mice with recurrent tumors were treated for a second round of Capecitabine when PDX reached a volume of over 200 mm3 (mice #35, #40 & #33). Mouse #40 did not respond to Capecitabine and the PDX specimen was extracted at 1080 mm3 and tagged as HBCx-95-CapaR. Samples studied in manuscript are indicated with a dashed rectangle. (b) Histograms of the distribution of scChIPseq raw and unique sequencing reads per cell in untreated HBCx-95 and Capecitabine-resistant HBCx-95-CapaR PDX. (c) Scatter plot displaying the log2 enrichments in cumulative single-cell versus bulk H3K27me3 ChIP-seq data, calculated based on the counts per million mapped reads in 50 kb genomic bins (n = 54,650 genomic bins). Pearson’s correlation score and P value are computed genome-wide. (d) Distribution of the log2 enrichments in 50 kb genomic bins for cumulative single-cell and bulk H3K27me3 ChIP-seq data. Amplitude between maximum and minimum log2 enrichments are indicated on the top right for each ChIP-seq method (either bulk or single-cell).

Supplementary Fig. 10 Detailed unsupervised analysis of scChIP-seq H3K27me3 profiles for mouse stromal cells from the HBCx-95 model.

(a) PCA plots with varying minimum number of unique mapped reads per cell showing that H3K27me3 scChIPseq principal components were mainly driven by cell coverage up to 1,600 reads per cell. (b) Hierarchical clustering and corresponding heatmap of cell-to-cell Pearson’s correlation scores based on their H3K27me3 profiles for n = 949 mouse stromal cells with a minimum coverage of 1,600 unique mapped reads. Cluster membership from hierarchical clustering, sample of origin (green for HBCx-95 and pink for HBCx-95-CapaR) and unique read count are indicated above the heatmap. (c) Distribution of the intra-cluster correlation scores for each cluster identified by hierarchical clustering. (d) Distribution of cell-to-cell Pearson’s correlation scores for real scChIP-seq datasets (in red) and randomized datasets (black curve). The upper percentile is indicated by a vertical black line. (e) Hierarchical clustering and corresponding heatmap of cell-to-cell Pearson’s correlation scores for ‘correlated’ cells only, representing 81.8% of all cells (see Methods). (f) Left: mean of all pairwise correlation scores between cluster’s members is plotted for k clusters ranging from 2 to 10. Right: mean intra-cluster correlation score for k clusters ranging from 2 to 10. At k = 3 clusters, the intra-cluster correlation is maximized. (g) Hierarchical clustering and corresponding heatmap of cell-to-cell consensus clustering score for scChIP-seq mouse stromal cells. Consensus score ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). Cluster membership is color coded above the heatmap.

Supplementary Fig. 11 Differential analysis of single-cell ChIP-seq H3K27me3 profiles of mouse stromal cells from the HBCx-95 model.

(a) Left panel: differentially bound regions in cells from each cluster versus all remaining cells identified by two-sided Wilcoxon signedrank test. Genomic regions were considered enriched (red) or depleted (green) if the adjusted P values were lower than 0.01 and the absolute log2 fold change greater than 1. Middle and right panels: barplots displaying the -log10 of adjusted P values from pathway analysis (using MSigDB, adjusted P value < 0.1) for regions depleted in H3K27me3 specifically in Chrom_c2 and Chrom_c3. The top 10 significant gene sets are indicated on the barplots. (b) Left panel: snapshot of a differentially depleted locus for cells from Chrom_c2 (Col4a1) versus all other cells, with cumulative single-cell profiles for each cell type and cluster (n = 276, 184 and 307 single-cells for Chrom_c1, Chrom_c2 and Chrom_c3 respectively). Differentially bound regions identified by two-sided Wilcoxon signed-rank test is indicated in gray with the corresponding adjusted P value and log2 fold change, calculated in respect to Chrom_c2. The number and the proportion of cells with H3K27me3 enrichment within each cluster are indicated above tracks. Right panels: t-SNE representation of scChIP-seq and scRNA-seq datasets, dots are colored according to H3K27me3 enrichment and expression signals in each cell for Col4a1 locus.

Supplementary Fig. 12 Clustering and differential analysis of single-cell RNA-seq profiles of mouse stromal cells from the HBCx-95 model.

(a) Hierarchical clustering and corresponding heatmap of cell-to-cell Pearson’s correlation scores. Cluster membership from consensus clustering, sample of origin (green for HBCx-95 [n = 617 single-cells] and pink for HBCx-95-CapaR [n = 926 single-cells]) and unique read count are indicated above the heatmap. (b) Consensus clustering scores for hierarchical clustering of stromal cells. Consensus score ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). Cluster membership is color coded above the dendrogram. (c) t-SNE plots with dots colored according to the expression signal of each marker gene, representative of each cell population, in each cell.

Supplementary Fig. 13 Clustering of single-cell ChIP-seq H3K27me3 profiles and scRNA-seq profiles of human tumor cells from the HBCx-95 model.

(a) Left: Plot of copy number in 0.5 Mb non-overlapping regions in Capecitabine-resistant PDX (HBCx-95-CapaR) versus untreated PDX (HBCx-95), obtained from the input of the bulk ChIP-seq experiments. Regions with a deviation greater than n = 2 standard deviations (dashed black lines) compared to the mean (diagonal black line) were removed for subsequent scChIP-seq analysis. Middle and right: snapshots of loci affected by copy number variation for bulk DNA profiles of Capecitabine-resistant PDX and untreated PDX indicated in gray. (b) Left panel: hierarchical clustering and corresponding heatmap of cell-to-cell Pearson’s correlation scores. Sample of origin (green for HBCx-95 [n = 474 single-cells] and pink for HBCx-95-CapaR [n = 838 single-cells]) and unique read count are indicated above the heatmap. Middle panel: t-SNE plots of scRNA-seq tumor cells, dots are colored according to the sample of origin and consensus clustering segmentation. Right panel: Consensus clustering scores for hierarchical clustering of scRNA-seq tumor cells. Consensus score ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). Cluster membership is color coded beneath the dendrogram. (c) Consensus clustering analysis for scChIP-seq dataset. Left panel: mean of all pairwise correlation score between cluster’s members is plotted for k clusters ranging from 2 to 10. At k = 2 clusters, the intra-cluster correlation is maximized. Right panel: hierarchical clustering and corresponding heatmap of cell-to-cell consensus clustering scores for scChIP-seq on tumor cells (HBCx-95 and HBCx-95-CapaR PDXs). Consensus scores ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). Cluster membership is color coded above heatmap. (d) Barplot displaying the -log10 of adjusted P values from pathway analysis for regions with depletion of H3K27me3 in resistant cells. The gene sets are indicated on the barplot. (e) Left panels: Aggregated H3K27me3 chromatin profiles for Chrom_c1 and Chrom_c2 (n = 457 and 794 single-cells for Chrom_c1 and Chrom_c2 respectively) are shown for the loci identified in Fig. 3f, as significantly differentially enriched and expressed (two-sided Wilcoxon signed-rank test). For each window indicated in gray, the log2 fold-change, the adjusted P value, the number and the proportion of cells with H3K27me3 enrichment within each cluster are indicated. Right panels: t-SNE representation of scRNA-seq datasets. Dots are colored according to expression signal in each cell.

Supplementary Fig. 14 Clustering of single-cell ChIP-seq profiles of human tumor cells from the HBCx-22 model.

(a) Histograms of the distribution of scChIP-seq raw and unique sequencing reads per cell in untreated HBCx-22 and Tamoxifen-resistant HBCx-22-TamR PDX. (b) Copy number in 0.5 Mb non-overlapping regions plotted for bulk DNA profiles of Tamoxifen-resistant PDX (HBCx-22-TamR) versus untreated PDX (HBCx-22). No aberrant variation in copy number was identified in this xenograft model between the resistant and sensitive tumors. The mean (diagonal black line) and deviations greater than n = 2 standard deviations compared to the mean (dashed black lines) are represented. (c) Consensus clustering analysis for scChIP-seq dataset. Left panel: mean of all pairwise correlation score between cluster’s members is plotted for k clusters ranging from 2 to 10. At k = 2 clusters, the intra-cluster correlation is maximized. Right panel: hierarchical clustering and corresponding heatmap of cell-to-cell consensus clustering scores for scChIP-seq on tumor cells (HBCx-22 and HBCx-22-TamR PDXs). Consensus scores ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). Cluster membership is color coded above heatmap. (d) Volcano plot representing adjusted P values (two-sided Wilcoxon signed-rank test) versus fold-changes for differential analysis comparing chromatin marks between Chrom_c2 and Chrom_c1 (n = 356 enriched and 137 depleted regions with adjusted P value < 0.01 and |log2FC| > 1). (e) Barplot displaying the -log10 of adjusted P values from pathway analysis for regions with depletion of H3K27me3 in cells from Chrom_c2. The gene sets are indicated on the barplot. (f) Hierarchical clustering and corresponding heatmap of cell-to-cell consensus clustering score for scRNA-seq tumor cells (HBCx-22 and HBCx-22-TamR PDXs). Consensus score ranges from 0 (white: never clustered together) to 1 (dark blue: always clustered together). Cluster membership is color coded above the heatmap. (g) Left panel: aggregated H3K27me3 chromatin profiles for Chrom_c1 and Chrom_c2 are shown for the ALCAM locus (n = 212 and 161 single-cells for Chrom_c1 and Chrom_c2 respectively). Differentially bound region identified by two-sided Wilcoxon signed-rank test is highlighted in gray with the corresponding adjusted P value and log2 fold-change. Middle panel: barplot displaying the proportion of cells with H3K27me3 enrichment in each cluster. The corresponding number of cells is indicated above the barplot. For each cluster, the origin of cells (green for HBCx-22 and pink for HBCx-22-TamR) is indicated below. Right panel: barplot displaying the average log2 fold-change for ALCAM expression level for cells in each cluster versus all remaining cells. The percentage of cells, within each cluster, with detectable ALCAM expression is indicated above the barplot. For each cluster, the origin of cells (green for HBCx-22 and pink for HBCx-22-TamR) is indicated below.

Supplementary information

Supplementary Information

Supplementary Figs. 1–14

Reporting Summary

Supplementary Tables

Supplementary Tables 1–5

Supplementary Video 1

Droplet re-injection and pairing.

Supplementary Video 2

Droplet fusion.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Grosselin, K., Durand, A., Marsolier, J. et al. High-throughput single-cell ChIP-seq identifies heterogeneity of chromatin states in breast cancer. Nat Genet 51, 1060–1066 (2019). https://doi.org/10.1038/s41588-019-0424-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41588-019-0424-9

This article is cited by

Search

Quick links

Nature Briefing: Cancer

Sign up for the Nature Briefing: Cancer newsletter — what matters in cancer research, free to your inbox weekly.

Get what matters in cancer research, free to your inbox weekly. Sign up for Nature Briefing: Cancer