Chromatin Velocity reveals epigenetic dynamics by single-cell profiling of heterochromatin and euchromatin

Abstract

Recent efforts have succeeded in surveying open chromatin at the single-cell level, but high-throughput, single-cell assessment of heterochromatin and its underlying genomic determinants remains challenging. We engineered a hybrid transposase including the chromodomain (CD) of the heterochromatin protein-1α (HP-1α), which is involved in heterochromatin assembly and maintenance through its binding to trimethylation of the lysine 9 on histone 3 (H3K9me3), and developed a single-cell method, single-cell genome and epigenome by transposases sequencing (scGET-seq), that, unlike single-cell assay for transposase-accessible chromatin with sequencing (scATAC-seq), comprehensively probes both open and closed chromatin and concomitantly records the underlying genomic sequences. We tested scGET-seq in cancer-derived organoids and human-derived xenograft (PDX) models and identified genetic events and plasticity-driven mechanisms contributing to cancer drug resistance. Next, building upon the differential enrichment of closed and open chromatin, we devised a method, Chromatin Velocity, that identifies the trajectories of epigenetic modifications at the single-cell level. Chromatin Velocity uncovered paths of epigenetic reorganization during stem cell reprogramming and identified key transcription factors driving these developmental processes. scGET-seq reveals the dynamics of genomic and epigenetic landscapes underlying any cellular processes.

Access options

Subscription info for Korean customers

We have a dedicated website for our Korean customers. Please go to natureasia.com to subscribe to this journal.

Rent or Buy article

Get time limited or full article access on ReadCube.

from$8.99

All prices are NET prices.

Data availability

Fastq files and raw count matrices have been deposited to the Array Express platform (https://www.ebi.ac.uk/arrayexpress/) with the following IDs: E-MTAB-9648, E-MTAB-10218, E-MTAB-2020, E-MTAB-10219, E-MTAB-9650, E-MTAB-9651 and E-MTAB-9659. Source data are provided with this paper.

Code availability

Code necessary to preprocess scGET-seq data is available at https://github.com/leomorelli/scGET (ref. 102) and https://github.com/dawe/scatACC (ref. 103). Illustrative code snippets for postprocessing are reported in Supplementary Data 2.

References

  1. 1.

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

    CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 

  2. 2.

    Greaves, M. Evolutionary determinants of cancer. Cancer Discov. 5, 806–821 (2015).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  3. 3.

    Liau, B. B. et al. Adaptive chromatin remodeling drives glioblastoma stem cell plasticity and drug tolerance. Cell Stem Cell 20, 233–246 (2017).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  4. 4.

    Hangauer, M. J. et al. Drug-tolerant persister cancer cells are vulnerable to GPX4 inhibition. Nature 551, 247–250 (2017).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  5. 5.

    Brock, A., Chang, H. & Huang, S. Non-genetic heterogeneity—a mutation-independent driving force for the somatic evolution of tumours. Nat. Rev. Genet. 10, 336–342 (2009).

    CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 

  6. 6.

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

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  7. 7.

    Sharma, S. V. et al. A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell 141, 69–80 (2010).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  8. 8.

    Flavahan, W. A., Gaskell, E. & Bernstein, B. E. Epigenetic plasticity and the hallmarks of cancer. Science 357, eaal2380 (2017).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  9. 9.

    Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat. Methods 10, 1213–1218 (2013).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  10. 10.

    Buenrostro, J. D. et al. Single-cell chromatin accessibility reveals principles of regulatory variation. Nature 523, 486–490 (2015).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  11. 11.

    Tatarakis, A., Behrouzi, R. & Moazed, D. Evolving models of heterochromatin: from foci to liquid droplets. Mol. Cell 67, 725–727 (2017).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  12. 12.

    Ninova, M., Tóth, K. F. & Aravin, A. A. The control of gene expression and cell identity by H3K9 trimethylation. Development 146, dev181180 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  13. 13.

    Nicetto, D. et al. H3K9me3-heterochromatin loss at protein-coding genes enables developmental lineage specification. Science 363, 294–297 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  14. 14.

    Nakayama, J., Rice, J. C., Strahl, B. D., Allis, C. D. & Grewal, S. I. Role of histone H3 lysine 9 methylation in epigenetic control of heterochromatin assembly. Science 292, 110–113 (2001).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  15. 15.

    Peters, A., O’Carroll, D. & Scherthan, H. Loss of the Suv39h histone methyltransferases impairs mammalian heterochromatin and genome stability. Cell 107, 323–337 (2001).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  16. 16.

    Patel, A. P. et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 1396–1401 (2014).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  17. 17.

    Aldridge, S. & Teichmann, S. A. Single cell transcriptomics comes of age. Nat. Commun. 11, 4307 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  18. 18.

    Henikoff, S., Henikoff, J., Kaya-Okur, H. & Ahmad, K. Efficient chromatin accessibility mapping in situ by nucleosome-tethered tagmentation. eLife 9, e63274 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  19. 19.

    Jacobs, S. A. & Khorasanizadeh, S. Structure of HP1 chromodomain bound to a lysine 9-methylated histone H3 tail. Science 295, 2080–2083 (2002).

    CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 

  20. 20.

    Lachner, M., O’Carroll, D., Rea, S., Mechtler, K. & Jenuwein, T. Methylation of histone H3 lysine 9 creates a binding site for HP1 proteins. Nature 410, 116–120 (2001).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  21. 21.

    Bannister, A. J. et al. Selective recognition of methylated lysine 9 on histone H3 by the HP1 chromo domain. Nature 410, 120–124 (2001).

    CAS 
    Article 

    Google Scholar 

  22. 22.

    Satpathy, A. T. et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nat. Biotechnol. 37, 925–936 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  23. 23.

    Cross, W. et al. The evolutionary landscape of colorectal tumorigenesis. Nat. Ecol. Evol. 2, 1661–1672 (2018).

    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  24. 24.

    Cross, W. et al. Stabilising selection causes grossly altered but stable karyotypes in metastatic colorectal cancer. Preprint at bioRxiv https://doi.org/10.1101/2020.03.26.007138 (2020).

  25. 25.

    Gézsi, A. et al. VariantMetaCaller: automated fusion of variant calling pipelines for quantitative, precision-based filtering. BMC Genomics 16, 875 (2015).

    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  26. 26.

    Misale, S. et al. Vertical suppression of the EGFR pathway prevents onset of resistance in colorectal cancers. Nat. Commun. 6, 8305 (2015).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  27. 27.

    Lupo, B. et al. Colorectal cancer residual disease at maximal response to EGFR blockade displays a druggable Paneth cell-like phenotype. Sci. Transl. Med. 12, eaax8313 (2020).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  28. 28.

    Laurent-Puig, P., Lievre, A. & Blons, H. Mutations and response to epidermal growth factor receptor Inhibitors. Clin. Cancer Res. 15, 1133–1139 (2009).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  29. 29.

    Wang, C. et al. Acquired resistance to EGFR TKIs mediated by TGFβ1/integrin β3 signaling in EGFR-mutant lung cancer. Mol. Cancer Ther. 18, 2357–2367 (2019).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  30. 30.

    Hu, T. & Li, C. Convergence between Wnt-β-catenin and EGFR signaling in cancer. Mol. Cancer 9, 236 (2010).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  31. 31.

    Sondka, Z. et al. The COSMIC Cancer Gene Census: describing genetic dysfunction across all human cancers. Nat. Rev. Cancer 18, 696–705 (2018).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  32. 32.

    Rondinelli, B. et al. Histone demethylase JARID1C inactivation triggers genomic instability in sporadic renal cancer. J. Clin. Invest. 125, 4625–4637 (2015).

    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  33. 33.

    Peric-Hupkes, D. et al. Molecular maps of the reorganization of genome–nuclear lamina interactions during differentiation. Mol. Cell 38, 603–613 (2010).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  34. 34.

    Hiratani, I. et al. Global reorganization of replication domains during embryonic stem cell differentiation. PLoS Biol. 6, 2220–2236 (2008).

    CAS 
    Article 

    Google Scholar 

  35. 35.

    Marchal, C. et al. Genome-wide analysis of replication timing by next-generation sequencing with E/L Repli-seq. Nat. Protoc. 13, 819–839 (2018).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  36. 36.

    Rondinelli, B. et al. H3K4me3 demethylation by the histone demethylase KDM5C/JARID1C promotes DNA replication origin firing. Nucleic Acids Res. 43, 2560–2574 (2015).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  37. 37.

    Wong, R. C. B. et al. L1TD1 is a marker for undifferentiated human embryonic stem cells. PLoS ONE 6, e19355 (2011).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  38. 38.

    Wong, Y. H. et al. Protogenin defines a transition stage during embryonic neurogenesis and prevents precocious neuronal differentiation. J. Neurosci. 30, 4428–4439 (2010).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  39. 39.

    Setty, M. et al. Characterization of cell fate probabilities in single-cell data with Palantir. Nat. Biotechnol. 37, 451–460 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  40. 40.

    Wang, C. et al. Reprogramming of H3K9me3-dependent heterochromatin during mammalian embryo development. Nat. Cell Biol. 20, 620–631 (2018).

    CAS 
    Article 

    Google Scholar 

  41. 41.

    Nicetto, D. & Zaret, K. S. Role of H3K9me3 heterochromatin in cell identity establishment and maintenance. Curr. Opin. Genet. Dev. 55, 1–10 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  42. 42.

    Burton, A. et al. Heterochromatin establishment during early mammalian development is regulated by pericentromeric RNA and characterized by non-repressive H3K9me3. Nat. Cell Biol. 22, 767–778 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  43. 43.

    Novo, C. L. et al. The pluripotency factor Nanog regulates pericentromeric heterochromatin organization in mouse embryonic stem cells. Genes Dev. 30, 1101–1115 (2016).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  44. 44.

    La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  45. 45.

    Wold, S., Sjöström, M. & Eriksson, L. PLS-regression: a basic tool of chemometrics. Chemom. Intell. Lab. Syst. 58, 109–130 (2001).

    CAS 
    Article 

    Google Scholar 

  46. 46.

    Eferl, R. et al. Development of pulmonary fibrosis through a pathway involving the transcription factor Fra-2/AP-1. Proc. Natl Acad. Sci. USA 105, 10525–10530 (2008).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  47. 47.

    Soares, E. & Zhou, H. Master regulatory role of p63 in epidermal development and disease. Cell. Mol. Life Sci. 75, 1179–1190 (2018).

    CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 

  48. 48.

    Zhu, M. & Zernicka-Goetz, M. Principles of self-organization of the mammalian embryo. Cell 183, 1467–1478 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  49. 49.

    Begley, C. G. et al. Molecular characterization of NSCL, a gene encoding a helix–loop–helix protein expressed in the developing nervous system. Proc. Natl Acad. Sci. USA 89, 38–42 (1992).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  50. 50.

    Lombardi, L. M. et al. MECP2 disorders: from the clinic to mice and back. J. Clin. Invest. 125, 2914–2923 (2015).

    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  51. 51.

    Martin Caballero, I., Hansen, J., Leaford, D., Pollard, S. & Hendrich, B. D. The methyl-CpG binding proteins Mecp2, Mbd2 and Kaiso are dispensable for mouse embryogenesis, but play a redundant function in neural differentiation. PLoS ONE 4, e4315 (2009).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  52. 52.

    Li, C. H. et al. MeCP2 links heterochromatin condensates and neurodevelopmental disease. Nature 586, 440–444 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  53. 53.

    Van Der Raadt, J., Van Gestel, S. H. C., Kasri, N. N. & Albers, C. A. ONECUT transcription factors induce neuronal characteristics and remodel chromatin accessibility. Nucleic Acids Res. 47, 5587–5602 (2019).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  54. 54.

    Rhee, H. S. et al. Expression of terminal effector genes in mammalian neurons is maintained by a dynamic relay of transient enhancers. Neuron 92, 1252–1265 (2016).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  55. 55.

    Cardoso-Moreira, M. et al. Gene expression across mammalian organ development. Nature 571, 505–509 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  56. 56.

    Kaya-Okur, H. S. et al. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nat. Commun. 10, 1930 (2019).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  57. 57.

    Wu, S. J. et al. Single-cell analysis of chromatin silencing programs in development and tumor progression. Preprint at bioRxiv https://doi.org/10.1101/2020.09.04.282418 (2020).

  58. 58.

    Stadhouders, R. et al. Transcription factors orchestrate dynamic interplay between genome topology and gene regulation during cell reprogramming. Nat. Genet. 50, 238–249 (2018).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  59. 59.

    Soufi, A., Donahue, G. & Zaret, K. S. Facilitators and impediments of the pluripotency reprogramming factors’ initial engagement with the genome. Cell 151, 994–1004 (2012).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  60. 60.

    Chen, J. Perspectives on somatic reprogramming: spotlighting epigenetic regulation and cellular heterogeneity. Curr. Opin. Genet. Dev. 64, 21–25 (2020).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  61. 61.

    Li, D. et al. Chromatin accessibility dynamics during iPSC reprogramming. Cell Stem Cell 21, 819–833 (2017).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  62. 62.

    Schwarz, B. A. et al. Prospective isolation of poised iPSC intermediates reveals principles of cellular reprogramming. Cell Stem Cell 23, 289–305 (2018).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  63. 63.

    Zviran, A. et al. Deterministic somatic cell reprogramming involves continuous transcriptional changes governed by Myc and epigenetic-driven modules. Cell Stem Cell 24, 328–341 (2019).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  64. 64.

    Lin, C., Ding, J. & Bar-Joseph, Z. Inferring TF activation order in time series scRNA-Seq studies. PLoS Comput. Biol. 16, e1007644 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  65. 65.

    Picelli, S. et al. Tn5 transposase and tagmentation procedures for massively scaled sequencing projects. Genome Res. 24, 2033–2040 (2014).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  66. 66.

    Machida, S. et al. Structural basis of heterochromatin formation by human HP1. Mol. Cell 69, 385–397 (2018).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  67. 67.

    Reznikoff, W. S. Transposon Tn5. Annu. Rev. Genet. 42, 269–286 (2008).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  68. 68.

    Zhu, Q. et al. BRCA1 tumour suppression occurs via heterochromatin-mediated silencing. Nature 477, 179–184 (2011).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  69. 69.

    Bertotti, A. et al. A molecularly annotated platform of patient-derived xenografts (‘xenopatients’) identifies HER2 as an effective therapeutic target in cetuximab-resistant colorectal cancer. Cancer Discov. 1, 508–523 (2011).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  70. 70.

    Reinhardt, P. et al. Derivation and expansion using only small molecules of human neural progenitors for neurodegenerative disease modeling. PLoS ONE 8, e59252 (2013).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  71. 71.

    Smith, T., Heger, A. & Sudbery, I. UMI-tools: modeling sequencing errors in unique molecular identifiers to improve quantification accuracy. Genome Res. 27, 491–499 (2017).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  72. 72.

    Lassmann, T. TagDust2: a generic method to extract reads from sequencing data. BMC Bioinformatics 16, 24 (2015).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  73. 73.

    Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Preprint at arXiv https://arxiv.org/abs/1303.3997 (2013).

  74. 74.

    Faust, G. G. & Hall, I. M. SAMBLASTER: fast duplicate marking and structural variant read extraction. Bioinformatics 30, 2503–2505 (2014).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  75. 75.

    Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A. & Manke, T. DeepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 42, 187–191 (2014).

    Article 
    CAS 

    Google Scholar 

  76. 76.

    Zhang, Y. et al. Model-based analysis of ChIP–seq (MACS). Genome Biol. 9, R137 (2008).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  77. 77.

    Breeze, C. E. et al. Atlas and developmental dynamics of mouse DNase I hypersensitive sites. Preprint at bioRxiv https://doi.org/10.1101/2020.06.26.172718 (2020).

  78. 78.

    Giansanti, V., Tang, M. & Cittaro, D. Fast analysis of scATAC-seq data using a predefined set of genomic regions. F1000Res. 9, 199 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  79. 79.

    Meuleman, W. et al. Index and biological spectrum of human DNase I hypersensitive sites. Nature 584, 244–251 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  80. 80.

    Quinlan, A. R. BEDTools: the Swiss-Army tool for genome feature analysis. Curr. Protoc. Bioinformatics https://doi.org/10.1002/0471250953.bi1112s47 (2014).

  81. 81.

    Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  82. 82.

    Polański, K. et al. BBKNN: fast batch alignment of single cell transcriptomes. Bioinformatics 36, 964–965 (2020).

    PubMed 
    PubMed Central 

    Google Scholar 

  83. 83.

    Traag, V. A., Waltman, L. & van Eck, N. J. From Louvain to Leiden: guaranteeing well-connected communities. Sci. Rep. 9, 5233 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  84. 84.

    Morelli, L., Giansanti, V. & Cittaro, D. Nested stochastic block models applied to the analysis of single cell data. Preprint at bioRxiv https://doi.org/10.1101/2020.06.28.176180 (2020).

  85. 85.

    Žitnik, M. & Zupan, B. Data fusion by matrix factorization. IEEE Trans. Pattern Anal. Mach. Intell. 37, 41–53 (2015).

    PubMed 
    Article 
    PubMed Central 

    Google Scholar 

  86. 86.

    Cho, S. W. et al. Promoter of lncRNA gene PVT1 is a tumor-suppressor DNA boundary element. Cell 173, 1398–1412 (2018).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  87. 87.

    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 (2009).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  88. 88.

    Zhao, H. et al. CrossMap: a versatile tool for coordinate conversion between genome assemblies. Bioinformatics 30, 1006–1007 (2014).

    PubMed 
    Article 
    CAS 
    PubMed Central 

    Google Scholar 

  89. 89.

    Karimzadeh, M., Ernst, C., Kundaje, A. & Hoffman, M. M. Umap and Bismap: quantifying genome and methylome mappability. Nucleic Acids Res. 46, e120 (2018).

    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 

  90. 90.

    Favero, F. et al. Sequenza: allele-specific copy number and mutation profiles from tumor sequencing data. Ann. Oncol. 26, 64–70 (2015).

    CAS 
    Article 

    Google Scholar 

  91. 91.

    Househam, J., Cross, W. C. H. & Caravagna, G. A fully automated approach for quality control of cancer mutations in the era of high-resolution whole genome sequencing. Preprint at bioRxiv https://doi.org/10.1101/2021.02.13.429885 (2021).

  92. 92.

    Caravagna, G., Sanguinetti, G., Graham, T. A. & Sottoriva, A. The MOBSTER R package for tumour subclonal deconvolution from bulk DNA whole-genome sequencing data. BMC Bioinformatics 21, 531 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  93. 93.

    Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. Preprint at arXiv https://arxiv.org/abs/1207.3907 (2012).

  94. 94.

    Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92 (2012).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  95. 95.

    Forbes, S. A. et al. COSMIC: mining complete cancer genomes in the catalogue of somatic mutations in cancer. Nucleic Acids Res. 39, 945–950 (2011).

    Article 
    CAS 

    Google Scholar 

  96. 96.

    Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  97. 97.

    Kaminow, B., Yunusov, D. & Dobin, A. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data. Preprint at bioRxiv https://doi.org/10.1101/2021.05.05.442755 (2021).

  98. 98.

    Harrow, J. et al. GENCODE: the reference human genome annotation for the ENCODE project. Genome Res. 22, 1760–1774 (2012).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  99. 99.

    Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281–291 (2019).

    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 

  100. 100.

    Kulakovskiy, I. V. et al. HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChIP–seq analysis. Nucleic Acids Res. 46, D252–D259 (2018).

    CAS 
    Article 

    Google Scholar 

  101. 101.

    Molineris, I., Grassi, E., Ala, U., Di Cunto, F. & Provero, P. Evolution of promoter affinity for transcription factors in the human lineage. Mol. Biol. Evol. 28, 2173–2183 (2011).

    CAS 
    PubMed 
    Article 

    Google Scholar 

  102. 102.

    Morelli, L. & Cittaro, D. scGET: pre-release of scGET repository. Zenodo https://doi.org/10.5281/zenodo.5095040 (2021).

  103. 103.

    Cittaro, D. scatACC: version 0.1. Zenodo https://doi.org/10.5281/zenodo.5095157 (2021).

Download references

Acknowledgements

We thank all the members of the COSR and Tonon laboratory for discussions, support and critical reading of the manuscript. We are grateful to E. Brambilla and F. Ruffini for preparation of the iPSCs and NPCs and A. Mira for assistance in the preparation of the organoids. We would like to thank S. de Pretis for the thoughtful discussions about chromatin velocity. We are grateful to G. Bucci for providing raw exome sequencing data and P. Dellabona for the coordination of the metastatic colon cancer sample collection and analysis. We also thank D. Gabellini, M. E. Bianchi, A. Agresti and S. Biffo for helpful discussions and for reviewing the manuscript. A.B. and L.T. are members of the EurOPDX Consortium. This work was partially supported by the Italian Ministry of Health with Ricerca Corrente and 5 × 1000 funds (S.M. and S.P.), by Associazione Italiana per la Ricerca sul Cancro (AIRC) investigator grants 20697 (to A.B.) and 22802 (to L.T.), AIRC 5 × 1000 grant 21091 (to A.B. and L.T.), AIRC/CRUK/FC AECC Accelerator Award 22795 (to L.T.), European Research Council Consolidator Grant 724748 BEAT (to A.B.), H2020 grant agreement 754923 COLOSSUS (to L.T.), H2020 INFRAIA grant agreement 731105 EDIReX (to A.B.), Fondazione Piemontese per la Ricerca sul Cancro-ONLUS, 5 × 1000 Ministero della Salute 2014, 2015 and 2016 (to L.T.), AIRC investigator grants (to G.T.) and by the Italian Ministry of Health with 5 × 1000 funds, Fiscal Year 2014 (to G.T.), AIRC 5 × 1000 ID 22737 (to G.T.) and the AIRC/CRUK/FC AECC Accelerator Award ‘Single Cell Cancer Evolution in the Clinic’ A26815 (AIRC number program 2279) (to G.T.).

Author information

Author notes

  1. Dalia Rosano

    Present address: Department of Surgery and Cancer, Imperial College London, London, UK

  2. These authors contributed equally: Martina Tedesco, Francesca Giannese.

Affiliations

  1. Università Vita-Salute San Raffaele, Milano, Italy

    Martina Tedesco, Paola Panina Bordignon & Gianvito Martino

  2. Functional Genomics of Cancer Unit, Division of Experimental Oncology, IRCCS San Raffaele Scientific Institute, Milano, Italy

    Martina Tedesco, Dalia Rosano, Oronza A. Botrugno & Giovanni Tonon

  3. Center for Omics Sciences, IRCCS San Raffaele Institute, Milano, Italy

    Francesca Giannese, Dejan Lazarević, Valentina Giansanti, Leonardo Morelli, Davide Cittaro & Giovanni Tonon

  4. Department of Informatics, Systems and Communication, University of Milano-Bicocca, Milano, Italy

    Valentina Giansanti

  5. Biochemistry and Structural Biology Unit, Department of Experimental Oncology, IEO, IRCCS European Institute of Oncology, Milano, Italy

    Silvia Monzani & Sebastiano Pasqualato

  6. Department of Oncology, University of Torino School of Medicine, Candiolo, Torino, Italy

    Irene Catalano, Elena Grassi, Andrea Bertotti & Livio Trusolino

  7. Candiolo Cancer Institute FPO- IRCCS, Candiolo, Torino, Italy

    Irene Catalano, Elena Grassi, Eugenia R. Zanella, Andrea Bertotti & Livio Trusolino

  8. Neuroimmunology Unit, Institute of Experimental Neurology, Division of Neuroscience, IRCCS San Raffaele Hospital, Milano, Italy

    Paola Panina Bordignon & Gianvito Martino

  9. Department of Mathematics and Geosciences, University of Trieste, Trieste, Italy

    Giulio Caravagna

  10. Hepatobiliary Surgery Division, IRCCS San Raffaele Hospital, Milano, Italy

    Luca Aldrighetti

Contributions

M.T. performed experiments and analyzed the data. F.G. devised the methodology and experimental design, performed experiments and analyzed data. D.L. devised the methodology. V.G. performed bioinformatic analysis. D.R. performed experiments and provided experimental assistance and expertise. L.M. performed bioinformatic analysis. S.M. performed cloning and transposase production. I.C. and E.R.Z. performed in vivo experiments. O.A.B. performed experiments related to culturing and maintenance of organoids. E.G. performed bioinformatic analysis. G.C. performed analysis on whole-exome data. P.P.B. designed and supervised the FIB reprogramming and iPSC differentiation experiments. A.B. designed and supervised in vivo experiments and reviewed the data. G.M. designed and supervised the FIB reprogramming and iPSC differentiation experiments and reviewed the data. L.A. provided the primary samples used for the organoid experiments. S.P. designed and supervised transposase production and reviewed the data. L.T. designed and supervised in vivo experiments and reviewed data. D.C. designed the study, performed bioinformatic analysis and wrote the manuscript. G.T. designed the study, analyzed data and wrote the manuscript.

Corresponding authors

Correspondence to
Davide Cittaro or Giovanni Tonon.

Ethics declarations

Competing interests

M.T., F.G., D.L., S.P., D.C. and G.T. have submitted a patent application, pending, covering TnH.

Additional information

Peer review information Nature Biotechnology thanks Kun Zhang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work

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

Extended data

Extended Data Fig. 1 Tn5 transposase is able to tagment compacted chromatin featuring H3K9me3.

a, General scheme of TAM-ChIP technique (created with BioRender.com). A primary antibody (ChIP-validated antibody, dark grey) binds to a specific histone modification (light grey) over the genome (blue-red). A secondary antibody (TAM-ChIP conjugate, blue) is linked to the Tn5 transposon, which is made of Tn5 transposase (yellow) and the respective barcoded adapters (green). Upon the binding of the secondary antibody to the primary antibody, the linked Tn5 transposase targets and cuts the genomic regions flanking the histone modification, adding the barcoded adapters. TAM-ChIP was performed on two biological replicates for each condition (H3K4me3, H3K9me3 and NoAb). b, H3K4me3 (green) and H3K9me3 (red) enrichment profiles obtained either by ChIP-seq or TAM-ChIP-seq, compared with Input ChIP control (violet). c, Enrichment profile of heterochromatic genes FAM5B, NTF3, CACNA1E obtained from TAM-ChIP libraries assessed by Real Time-qPCR confirms Tn5 is able to target heterochromatic loci when redirected by H3K9me3 antibody. For each biological replicate three technical replicates were analyzed by Real-Time qPCR; one of the two H3K4me3 biological replicate was excluded because no appreciable signal was detected for any condition. Whiskers represent standard deviations (n = 3 technical replicates). Data shown in b and c refer to experiments performed on Caki-1 cell line.

Source data

Extended Data Fig. 2 Hybrid CD (HP1α)-Tn5 targets H3K9me3 chromatin regions.

a, Two different lengths of HP1α polypeptide (spanning amino acids 1-93 and 1-112) were linked to Tn5, using either a 3 or 5 poly-tyrosine–glycine–serine (TGS) linker, resulting in four hybrid construct, TnH#1-4. TnH#1 made of 1-93aa (HP1α) – 3x(TGS) – Tn5; TnH#2 made of 1-93aa (HP1α) – 5x(TGS) – Tn5; TnH#3 made of 1-112aa (HP1α) – 3x(TGS) – Tn5; TnH#4 made of 1-112aa (HP1α) – 5x(TGS) – Tn5. The 1-93 or 1-112aa spanning regions of HP1α include 1-75aa of CD followed by 18 or 37aa of natural linker, respectively (Created with BioRender.com). b-c, Tagmentation profiles relative to the four hybrid constructs (TnH#1-4) showing no difference in tagmentation efficiency relative to the native Tn5 enzyme (Nextera and Tn5 in-house produced) when targeting either genomic DNA, panel b, or native chromatin on permeabilized nuclei, panel c. d, Enrichment profiles relative to ATAC-seq performed with the four hybrid constructs (TnH#1-4, red) compared with native Tn5 enzyme (Nextera and Tn5 in-house produced) and with H3K4me3 and H3K9me3 ChIP-seq signals (green). e, Distribution of the enrichment of four TnH hybrid constructs (TnH#1-4) relative to genomic background in regions enriched for H3K4me3 (orange) or H3K9me3 (blue) expressed as log2(ratio) of the signal over the genomic Input. Enrichment over the same regions for native Tn5 enzyme (Nextera and Tn5 in-house produced), H3K4me3 and H3K9me3 ChIP-seq are reported as reference. Ec: global enrichment over H3K9me3-marked regions; Eo: global enrichment over H3K4me3-marked regions; Mc: modal enrichment over H3K9me3-marked regions; Mo: modal enrichment over H3K4me3-marked regions. Data shown in b, c and d refer to experiments performed on Caki-1 cell line.

Extended Data Fig. 3 Optimization of ATAC-seq protocol introducing a combination of Tn5 and TnH transposases.

a, Effect of altering Tn5 (green) to TnH (red) ratio on tagmentation profiles when adding both enzymes simultaneously at the beginning of the 60 minutes of the transposition reaction. b, Sequential addition of the same quantity of Tn5 and then TnH enzyme after 30 minutes of the transposition reaction results in a balanced distribution of enrichment signals between the two enzymes. Experiments performed on Caki-1 cell line.

Extended Data Fig. 4 Characteristic of scGET-seq data.

a Abundance of unique cell barcodes retrieved by scATAC-seq performed on Caki-1 cells using the provided ATAC transposition enzyme (10X Tn5; 10X Genomics) (blue) compared to cell barcodes countable by TnH (orange) or Tn5 (green) alone. scGET-seq performance (Tn5 + TnH) is represented in red. The curves are largely overlapping, indicating no evident bias in single cell identification; b Distribution of per-cell normalized coverage over fixed-size genomic bins (5 kb) is reported for 10X Tn5 (blue) and for signal obtained by TnH (orange) and Tn5 (green). While Tn5 is comparable to 10X Tn5, TnH returns higher and less overdispersed per-bin coverages. White dot in boxplots reprents the median, boxes span between the 25th and 75th percentiles, whiskers extend 1.5 times the interquartile range. n = 3363, 1281 and 1537 cells in one experiment; c Saturation analysis for selected libraries. Dotted lines show the fitted incomplete Gamma functions on subsampled data; red solid lines show subsampling data from the same libraries; d Tn5 (green) and TnH (red) enrichment profiles obtained from scGET-seq (pseudo-bulk) or from ATAC-seq performed by using the two enzymes separately, compared with H3K4me3 (green) and H3K9me3 (red) ChIP-seq data. Data shown refer to experiments performed on Caki-1 cells.

Extended Data Fig. 5 Copy Number analysis at multiple resolutions.

a, Segmentation profiles in individual cells profiled by 10X Tn5 (scATAC-seq) (left panel) or TnH scGET-seq (right panel) at 500 kb. b, Segmentation profiles in individual cells profiled by 10X Tn5 (scATAC-seq) (left panel) or TnH scGET-seq (right panel) at 1 Mb. c, Segmentation profiles in individual cells profiled by 10X Tn5 (scATAC-seq) (left panel) or TnH scGET-seq (right panel) at 10 Mb. On top of each heatmap the genome-wide coverage of bulk sequencing of corresponding cell lines is represented. Centromeric regions and gaps (in white) have been excluded from the analysis.

Extended Data Fig. 6 Characterization of Patient Derived Organoids.

a, evaluation of clonal structure of two PDO (CRC6 and CRC17) by exome sequencing; the histogram show the distribution of the cancer cell fraction estimated from the analysis of somatic mutations; in both organoids we observe a monoclonal structure b, 5X (left panel) and 10X (right panel) magnification contrast phase images of PDO #CRC17 obtained from a liver metastasis of a CRC patient (n>5); c absolute copy number of CRC17 and CRC6 as revealed by whole exome sequencing; data in panel c are equivalent to barplots over heatmaps in Fig. 3a.

Extended Data Fig. 7 scGET-seq analysis on PDX samples.

a, UMAP embedding of individual cells as in Fig. 3, colored by the time PDX were harvested. b, Segmentation profiles in individual cells profiled by scGET-seq at 1 Mb resolution expressed as log2(ratio) over the median signal. Cells are clustered according to genetic clones. Red: positive values; Blue: negative values. Centromeric regions (white) have been excluded from the analysis because they correspond to low mapping and not fully characterized regions.

Extended Data Fig. 8 scGET-seq profiling of NIH-3T3 cells knocked-down for Kdm5c.

a, Distribution of early-to-late ratio of 2-stage Repli-seq data for NIH-3T3 cells. Violin plots represent the value of log2(E/L) values over DHS regions which are differential in the high-vs-low coverage cells in Fig. 4a (Mann-Whitney U= 36169.5, p = 1.403e-84). White dot in boxplots represents the median, boxes span between the 25th and 75th percentiles, whiskers extend 1.5 times the interquartile range. n = 35438 regions. b, Distribution of lamin-B1 DamID scores for NIH-3T3 cells. Violin plots represent the value of DamID scores over DHS regions which are differential in the high-vs-low coverage cells in Fig. 4a (Mann-Whitney U = 723.0, p = 4.621e-6). White dot in boxplots represents the median, boxes span between the 25th and 75th percentiles, whiskers extend 1.5 times the interquartile range. n = 35438 regions. c, UMAP embedding of individual cells coloured by cell groups, identified by Leiden algorithm with resolution parameter set to 0.2. d, Results of the linear model calculating the group-wise differences between TnH and Tn5 enrichment. For each group we reported the coefficient of the model, the p-value and the Benjamini-Hochberg corrected p-value. Values are reported for the two genomic regions including the Major primers (see text). Barplot indicates the proportion of shScr-treated for each cell group.

Extended Data Fig. 9 scGET-seq profiling of a developmental model of iPSC.

a, UMAP embedding of individual cells colored by the probability of being included in a trajectory branch estimated by Palantir. Three major branches have been identified, roughly corresponding to the three cell types profiled in this study. b, Schematic representation of the phase portraits underlying Chromatin Velocity. In RNA-velocity, the time derivative of the unspliced/spliced RNA is used to estimate synthesis or degradation of RNA; in Chromatin Velocity, the same procedure is applied on Tn5/TnH data to estimate chromatin relaxation or compaction. d, UMAP embedding of individual cells colored by cell clusters. e, Heatmap shows average expression profiles of TF with the top 10 most negative on PLS2 during the early brain development. Darker color indicates higher expression. w.p.c.: weeks post conception.

Supplementary information

Supplementary Table 1

Counts of cells from organoid CRC6 or CRC7 found in different clones identified using TnH (above) or Tn5 (below).

Supplementary Table 2

Enrichment analysis over KEGG pathways and Reactome pathways of genes associated with DHS sites that are found to be differentially enriched in epigenetic clones. Enrichment was performed using the Enrichr platform.

Supplementary Table 3

Mutations: list of somatic mutations of the primary tumor as a result of exome sequencing data. scGET-seq mutations: list of mutations profiled by scGET-seq. Only variants that have an impact on protein sequence have been reported.

Supplementary Table 4

Analysis of differential Tn5 signal enrichment according to different cell types. For each cell type, we report log fold change, P value and adjusted P value as a result of a t-test over each region. For each region, we report the closest genes (GENCODE v36) and the distance. We also report the log fold change, P value and adjusted P value of differential expression of the associated genes in each cell type

Supplementary Table 5

Analysis of differential Tn5 signal enrichment with respect to the cell entropy as estimated by Palantir. Regions are sorted for decreasing coefficient of the generalized linear model. Genes associated with regions by proximity are also reported.

Supplementary Table 6

Enrichment analysis of genes associated with top DHS regions with the dynamical profile. Analysis was performed using gProfiler.

Supplementary Table 7

Analysis of global transcription factor activity. HOCOMOCO v11 ID, PWM identification code; Gene Symbol, associated gene symbol; PLS1 and PLS2, loading of the TF after PLS regression, corresponds to the horizontal/vertical displacement of the TF arrows in Fig. 6e.

Supplementary Table 8

Sequencing statistics for all scGET-seq experiments presented in the manuscript. n_reads, number of sequencing fragments; n_reads_in_cell, number of fragments associated to a cell; n_duplicated, number of PCR duplicates; target cells, number of target cells in the experiment; PF cells, number of cells passing the initial processing filters (coverage by cell and by region); Compound Coverage, coverage estimate as number of mapped reads in cells (without duplicates) by read length divided by genome size; Per cell Coverage, average per cell coverage as Compound Coverage divided by the number of PF cells.

Supplementary Data 1

Amino acid sequences of TnH constructs (TGS residues underlined; H stands for histidine residue that is an artifact introduced as a consequence of the cloning strategy); Modified Tn5ME-A and TnHMe-A sequences with Tn- or TnH-associated barcode are underlined.

Supplementary Data 2

Representative code snippets to postprocess scGET-seq data.

About this article

Verify currency and authenticity via CrossMark

Cite this article

Tedesco, M., Giannese, F., Lazarević, D. et al. Chromatin Velocity reveals epigenetic dynamics by single-cell profiling of heterochromatin and euchromatin.
Nat Biotechnol (2021). https://doi.org/10.1038/s41587-021-01031-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1038/s41587-021-01031-1

Read More

Written by 

Leave a Reply

Your email address will not be published. Required fields are marked *