Single-cell EpiChem jointly measures drug–chromatin binding and multimodal epigenome – Nature Methods

Oligonucleotides

The sequences of oligonucleotide primers used in this work are available in Supplementary Table 1.

Chemical synthesis

The synthesis of Dox-btn was based on the previously reported method by the BioDuro-Sundia company8. NMR spectra were recorded on a Bruker 400 MHz Advance III Spectrometer for 1H NMR in DMSO-d6 and analyzed using MestReNova 14.2.3 software. NMR data were reported as follows: chemical shifts in parts per million (ppm) referring to the solvent residual peak, multiplicities (s, singlet; d, doublet; t, triplet; q, quartet; m, multiplet; br, broad) and coupling constant (J) values in Hz. LC–MS was performed on an Agilent ESI-MS (Agilent) connected to an Agilent 1260 Infinity system. Dox (T1020, CAS: 25316-40-9), JQ1 (T2110, CAS: 1268524-70-4) and THZ1 (T3664, CAS: 1604810-83-4) were obtained from TargetMol. JQ1-btn (HY-145667, CAS: 1635437-52-3) and THZ1-btn (HY-128867, CAS: 1604811-14-4) were obtained from MedChemExpress53,54.

Human participants

This study was approved by the Research Ethics Committee of the Institute of Basic Medical Sciences, Chinese Academy of Medical Sciences, School of Basic Medicine Peking Union Medical College (2022173). Written informed consent was obtained from each patient. Fresh tumor tissues were collected from treatment-naive patients with CRC who underwent primary curative resection within 30 min after the operation.

Cell and organoids culture

The mES cells (original from Novus Biologicals, NBP1-41162) were provided by H. Deng, Peking University, and maintained in this laboratory. The HGC27 cell line (original from ECACC, CVCL_1279) and SNU16 cell line (original from ATCC, CRL-5974) were provided by S. Shu, Peking University Cancer Hospital & Institute, and maintained in this laboratory. The human K562 cell line (original from ATCC, CCL-243) was maintained in this laboratory. All cells were cultured at 37 °C with 5% CO2. K562 cells were cultured in RPMI 1640 (Gibco) supplemented with 10% fetal bovine serum (FBS; Sigma). SNU16 cells were cultured in DMEM (Hyclone) supplemented with 10% FBS. HGC27 cells were cultured in RPMI 1640 (Gibco) supplemented with 20% FBS. The mES cells were cultured in DMEM (Hyclone) supplemented with 15% FBS, 1% Glutamax (Gibco), 0.1 mM 2-mercaptoethanol (Sigma), 1% MEM nonessential amino acids (Cellgro), 1% nucleoside (Millipore) and 1,000 U ml−1 recombinant leukemia inhibitory factor (Millipore). CRC organoids were established as described previously and cultured in CRC organoid medium55.

Cell imaging

SNU16, K562 and HGC27 cells were treated with JQ1 (2 µM), JQ1-btn (2 µM), THZ1 (1 µM), THZ1-btn (1 µM), Dox (3 µM) or Dox-btn (3 µM) for 3 h, respectively. Cells were fixed with 0.25% formaldehyde for 5 min on ice. CRC organoids were treated with JQ1 (5 µM), JQ1-btn (5 µM), THZ1 (5 µM), THZ1-btn (5 µM), Dox (3 µM) or Dox-btn (3 µM) for 6 h. Subsequently, the organoids were fixed with 4% formaldehyde for 10 min. After washing twice with 0.1% BSA/PBS, cells were dropped onto the slides for 2 h at room temperature. Slides were then permeabilized with 0.5% TX-100 in PBS and blocked with 3% BSA/PBS for 30 min. Finally, slides were incubated overnight at 4 °C with primary rabbit monoclonal anti-biotin (CST, 5597S, 1:250 dilution), and goat anti-rabbit Alexa Fluor 488 secondary antibodies (4412S, CST) were used for visualization of drug. Cell nuclei were stained with 4,6-diamidino-2-phenylindole (C1002, Beyotime), and the slides were scanned under a confocal microscope.

Drug treatment

CRC organoids were cultured in organoid medium supplemented with 5% Matrigel and plated in low-attachment six-well plates. After a 24-h recovery period, the organoids were exposed to either DMSO or small-molecule drugs at their respective IC50 concentrations for 3 or 5 days. Following the treatment, the organoids treated with DMSO or small-molecule drugs were collected and subjected to single-cell EpiChem assay for further analysis.

Cell viability assays

For dose–response assays, drugs were serially diluted in medium. These diluted compounds were then added to 96-well plates, which were initially seeded with 1.5 × 103 cells per well, resulting in a final concentration of 0.1% DMSO. After incubating drug for 96 h, CellTiter-Glo Luminescent Cell Viability Assay (G7572, Promega) was added to each well. The plates were further incubated for 10 min, and luminescence was measured using a SYNERGY H1 microplate reader (BioTek). CRC organoids suspended in 5% Matrigel containing organoid culture medium were plated in 384-well plates at a density of 1,000 per well. Gradient dilutions of drugs were added to the respective wells in the 384-well plates. After 5 days of treatment, the plates were equilibrated to room temperature and cell viability was assessed using the CellTiter-Glo viability assay (Promega), following the instructions provided by the manufacturer.

Barcoded protein A-Tn5 preparation

The production of PAT was conducted in-house according to the previous protocol12. In brief, pET28a-His-pA-Tn5 was transformed into BL21 (DE3)-competent cells. The colonies obtained were then cultured in 500 ml LB medium. To induce PAT expression, 0.2 mM IPTG was added when the optical density of the bacterial culture reached 0.5–0.8. Incubation of the culture was carried out at 23 °C and 80 rpm for 5 h. The bacterial pellets were collected and lysed using HGX buffer, followed by sonication. Genomic DNA in the lysate was precipitated with polyethyleneimine (Sigma). The clear lysate was then applied to a nickel-nitrilotriacetic acid column (QIAGEN) and PAT was eluted from the column using 100 mM imidazole. The concentration of PAT was quantified by Coomassie blue staining after resolution in a 7.5% PAGE gel. Next, PAT was mixed with the annealed barcoded adaptor (Supplementary Table 1) at an equal molar ratio of 37.5 µM, and the mixture was incubated at 25 °C for 1 h. The assembled PAT was either stored at −20 °C or used directly for further experiments.

Sample fixation

Cell pellets were resuspended in 1 ml 0.1% BSA/PBS, before 7 μl 36.5% formaldehyde was added for incubation on ice for 5 min. To terminate fixation, 14 μl 2.5 M glycine was added and the tubes were inverted several times and kept on ice for 5 min. Subsequently, cells were collected at 300g for 3 min at 4 °C and washed with 1 ml 0.1% BSA/PBS twice. Finally, cells were resuspended in cold 1% BSA/PBS and nine volumes of cold methanol were added dropwise on ice for storage at −80 °C.

EpiChem

In vivo

Cells were treated with 10 µM JQ1-btn for 15 min first; then cells were collected and fixed. Chromatin interactions of small molecules were captured by anti-biotin antibodies in antibody buffer (20 mM HEPES, pH 7.5, 150 mM NaCl, 0.5 mM Spermidine (Sigma), 1× Protease Inhibitor Cocktail (Roche), 10 mM sodium butyrate, 0.5 mM EDTA, 0.01% Triton X-100 and 0.01% digitonin) at 4 °C overnight. Then targeted tagmentation and library preparation were undertaken as follows.

In vitro

Cells after fixation were incubated with 2.5 µM JQ1-btn with anti-biotin antibodies in antibody buffer at 25 °C for 1 h. After washing in DIG wash buffer twice (20 mM HEPES, pH 7.5, 150 mM NaCl, 0.5 mM Spermidine, 10 mM sodium butyrate and 0.01% digitonin), cells were resuspended in 100 µl antibody buffer with secondary antibodies (Invitrogen, A32790, 1:500 dilution), followed by incubation at 4 °C for 30 min. Cells were washed twice in DIG wash buffer and incubated 9 µg ml−1 PAT–MEA and 9 µg ml−1 PAT–MEB in DIG-300 wash buffer (20 mM HEPES, pH 7.5, 300 mM NaCl, 0.5 mM Spermidine, 1× Protease Inhibitor Cocktail, 10 mM sodium butyrate, 0.01% Triton X-100 and 0.01% digitonin). Cells were washed twice in DIG-300 wash buffer and resuspended in 50 µl reaction buffer. The reaction was initiated at 37 °C for 1 h and terminated with DIG-300 wash buffer with 5 mM EDTA at room temperature for 5 min. Then, 2,000 cells were washed and resuspended in 4 µl Lysis Buffer (10 mM Tris-HCl, pH 8.5, 0.05 % SDS and 0.1 mg ml−1 Proteinase K). Following lysis at 55 °C for 15 min, 1 µl 10 mM phenylmethylsulfonyl fluoride (PMSF) and 1 µl 1.8% Triton X-100 were added to each well and the plate was incubated at 37 °C for 10 min to quench SDS. A total of 50 µl of PCR mix was added to the same tube comprising 25 µl 2× Tag Master mix, 0.5 µl 10 µM Nextera i7 primer, 0.5 µl 10 µM Nextera i5 primer, 1 µl 25 mM MgCl2 and 16 µl ddH2O. PCR enrichment was then conducted with one cycle of 72 °C for 5 min, one cycle of 95 °C for 2 min and 11 cycles of 98 °C for 30 s, 63 °C for 30 s and 72 °C for 1 min, followed by one cycle of 72 °C for 7 min and a hold at 4 °C. Afterward, 50 µl (0.9×) of custom AMPure XP beads were added to each well, mixed thoroughly, and DNA was purified and eluted using 10 µl ddH2O. For double-size selection purification, 0.5× + 0.5× of AMPure beads were used. DNA was purified and eluted with 20 µl ddH2O.

Single-cell EpiChem

Probe–primary antibody–PAT T7 complex incubation and tagmentation

A small-molecule primary Ab–PAT T7 complex was assembled as previously described10,28 with a few modifications. Compound stock solution in DMSO was diluted to 10 µM in probe solution (2 mM EDTA, 0.1% BSA (Sigma) and 0.05% digitonin in wash buffer). Then, 3.34 μl probe solution (10 µM), 0.5 μg antibody (3.33 pmol), 0.22 μl pre-assembled T7-barcoded PAT (8.25 pmol) and 5 μl antibody buffer were mixed thoroughly and incubated at 25 °C for 1 h. Cells were resuspended in 95 µl antibody buffer with the addition of 5 µl of primary Ab–PAT T7 complex and incubated at 4 °C for 4 h (For JQ1-btn, the final concentration was diluted to 2.5 µM and incubated at 25 °C for 1 h. For THZ1-btn, the final concentration was diluted to 0.5 µM and incubated at 4 °C for 4 h. For Dox-btn, the final concentration was diluted to 0.25 µM and incubated at 25 °C for 1 h). Cells were washed three times with DIG-300 wash buffer, and were resuspended in 50 µl reaction buffer consisting of 10 mM TAPS-NaOH (pH 8.3), 5 mM MgCl2, 1 mM dithiothreitol, proteinase inhibitor and 10 mM sodium butyrate. The reaction was initiated at 37 °C for 1 h. The reaction was stopped by incubating cells in 180 μl DIG-300 wash buffer with 5 mM EDTA for 5 min.

Primary antibody–PAT T7 complex incubation and tagmentation

The 0.5 μg antibody (3.33 pmol), 0.22 μl pre-assembled T7 barcoded PAT (8.25 pmol) and 5 μl wash buffer were mixed thoroughly and incubated at 25 °C for 1 h. After washing three times with DIG-300 wash buffer, cells were resuspended in 95 µl antibody buffer with 5 μl primary Ab–PAT T7 complex and incubated at 25 °C for 1 h. Cells were washed three times with DIG-300 wash buffer and resuspended in 50 µl reaction buffer. The reaction was initiated at 37 °C for 1 h. The reaction was stopped by incubating cells in 180 μl DIG-300 wash buffer with 5 mM EDTA at room temperature for 5 min.

Secondary antibody incubation and tagmentation

A secondary antibody diluted at 1:500 was added to 100 µl antibody buffer, followed by incubation at 4 °C for 15 min. After washing twice with DIG-300 wash buffer, cells were resuspended in 100 µl antibody buffer with 9 µg ml−1 of barcoded PAT-T5, followed by incubation at 4 °C for 30 min. Cells were washed twice with DIG-300 wash buffer and resuspended in 50 µl reaction buffer. The reaction was initiated at 37 °C for 1 h and terminated with 180 μl DIG-300 wash buffer with 5 mM EDTA for 5 min.

Chromatin tagmentation

Chromatin tagmentation was carried out as previously described56. Cells were suspended in 50 μl reaction buffer with 2.5 μM Tn5–T5/T7 transposome complex. The mixture was then incubated at 30 °C for 30 min and terminated with 180 μl DIG-300 wash buffer with 5 mM EDTA for 5 min.

Hybridizations and ligation

Ligation-based barcoding was carried out as previously described with a few modifications27. Cells were washed twice with NSB buffer and resuspended in 4.5 ml Hybridization Mix (1× T4 ligation buffer, 0.05% Triton X-100 and 0.25× NSB). Then, 40 µl cells were added to each of the 96 wells in the first-round barcoding plate, which already contained 10 µl annealed DNA barcodes and was incubated for 30 min at 25 °C (300 rpm). Then, 10 µl of round 1 blocking oligonucleotides was added into the plate and incubated for 30 min at 25 °C (300 rpm). Cells from all 96 wells were combined and 50 µl cells were added to each of the 96 wells in the second-round barcoding plate and incubated for 30 min at 25 °C (300 rpm). Then, 10 µl of round 2 blocking oligonucleotides was added into the plate, and incubated for 30 min at 25 °C (300 rpm). Cells from all wells were combined and washed twice with NSB buffer and centrifuged at 1,000g for 3 min. Subsequently, cells were resuspended in 200 µl ligation mix composed of 1× T4 ligation buffer, 20 U µl−1 T4 DNA ligase, 0.05% Triton X-100 and 0.2× NSB, and incubated for 30 min at 25 °C (300 rpm).

Redistributing cells and releasing DNA

Cells were washed and filtered through a 70-µm cell strainer to remove cell clumps. Subsequently, 1,000–3,000 cells were added into each well of new 96-well plates that contained 4 µl of Lysis Buffer. Following brief lysis at 55 °C for 15 min, 1 µl 10 mM PMSF and 1 µl 1.8% Triton X-100 were added to each well and the plate was incubated at 37 °C for 10 min to quench SDS.

Library preparation

A total of 50 µl of PCR mix was added to each well, comprising 10 µl 5× KAPA HiFi buffer, 1 µl 10 mM dNTP Mix, 2 µl 10 µM P7 connector primer, 2 µl 10 µM TruSeq P5 primer, 0.5 µl 1 U µl−1 KAPA HiFi HotStart DNA polymerase, 1 µl 25 mM MgCl2 and 26.5 µl ddH2O. PCR enrichment was then conducted with one cycle of 72 °C for 5 min, one cycle of 95 °C for 3 min, and 12 cycles of 98 °C for 20 s, 65 °C for 30 s and 72 °C for 30 s, followed by one cycle of 72 °C for 5 min and a hold at 4 °C. Afterward, 0.9× custom AMPure XP beads were used for DNA purification and followed by the addition of 50 µl PCR mix containing 10 µl 5× KAPA HiFi, 1 µl 10 mM dNTP Mix, 2.5 µl 10 µM P7 primer, 2.5 µl 10 µM P5 primer and 0.5 µl 1 U µl−1 KAPA HiFi HotStart DNA polymerase. The second PCR enrichment was carried out with one cycle at 72 °C for 5 min, one cycle at 98 °C for 3 min, six cycles at 98 °C for 20 s, 65 °C for 30 s, 72 °C for 1 min, and one cycle at 72 °C for 1 min, followed by hold at 4 °C. Finally, the library was purified once with 0.5× + 0.4× AMPure XP beads.

Data processing

scEpiChem data processing

scEpiChem data were processed to generate unique and non-duplicated reads as previously described12. In brief, we evaluated the quality of scEpiChem sequencing data by FastQC (v.0.11.5). Low-quality bases and adaptors were removed by Cutadapt (v.1.11) with the following parameters: -q 20 -O 10–trim-n -m 30–max-n 0.1. Clean paired-end EpiChem reads were then mapped to the human and mouse reference genome hg19 and mm10 using Bowtie2 (v.2.2.9)57. The mapped reads with MAPQ greater than 30 were considered as uniquely mapped reads, which were sorted using SAMtools (v.1.9) and used for subsequent analyses. PCR duplicates were removed by Picard (v.2.2.4) (http://broadinstitute.github.io/picard) with default parameters. Only uniquely mapped, non-duplicated reads were used for peak calling by MACS2 (v.2.1.1)58 and the following analyses.

Visualization and correlation analysis of scEpiChem data

We used the bamCoverage function in deepTools (v.2.2.3)59 to calculate genome coverage in BAM files and generated track files (bigWig format) for scEpiChem data. For the visualization of the signals of histone modifications/drug binding at specific loci, bigwig files were uploaded into Integrative Genomics Viewer (v.2.3)60 together with corresponding reference tracks. To evaluate correlations of JQ1-btn and BRD4 in K562 and HGC27, we calculated the normalized average scores in Chem-map BRD4 peak regions (51,963 peaks)10. The correlation was calculated between different groups and plotted by deepTools (v.2.2.3).

Collision rate

To evaluate the efficiency of split-and-pool strategy for single-cell labeling of scEpiChem, we generated scatter-plots using the proportion of reads mapped to human or mouse genome in each barcode combination in Fig. 1. Barcodes with <95% of aligned reads mapped to one species were classified as collisions61.

Single-cell scEpiChem data processing

The scEpiChem data were demultiplexed by custom scripts. In brief, the analysis pipeline of scEpiChem data processing consisted of following steps: (1) creating the single-cell whitelist using UMI-tools (v.1.1.2)62 based on the molecular design of scEpiChem with the parameter ‘–set-cell-number=n’; (2) extracting the paired-end reads based on the single-cell whitelist; (3) mapping to the human or mouse reference genomes by bowtie2 (v.2.2.9); (4) keeping unique mapped reads and removing PCR duplicates; and (5) adding cell barcode information to the BAM files and generating single-cell BAM files63.

Single-cell quality control and filtering

After removing doublets or multiplets in each batch by DoubletFinder (doubletFinder_v.3)64, we retained ~85% single cells (>400 unique fragments) with tri-omics information subjected to downstream analysis.

FRiP score

For each single-cell BAM file, the FRiP score calculation proceeds as follows: conversion of BAM to BED format using bedtools bamtobed, followed by intersection of reads with the peak regions using bedtools intersect. This tool is employed to count the number of reads that overlap with any of the defined peak regions, using the BED format reads and the peak file as input. Calculation of total read count within each BAM file was performed using SAMtools view -c. The FRiP score for each cell was then calculated as the ratio of reads intersecting with peaks to the total number of reads, expressed as a decimal to four significant figures. The computed FRiP scores for all analyzed cells were collected and written to a text file, preserving both the cell identifier and the corresponding FRiP score.

Dimensionality reduction and visualization of tri-omics in single cells

Single-cell data were processed using cisTopic v.3 (ref. 65) to generate a cell-peak matrix, and further processed using Signac (v.1.9.0)66 to perform term frequency-inverse document frequency (TF-IDF) normalization and latent semantic indexing clustering. The gene ACTIVITY assay in Seurat Object of cell-peak or cell-bin matrix for each small molecular and histone modification was generated by GeneActivity function of Seurat v.4 (ref. 38) with default parameters. The Seurat Object of AVTIVITY assay was processed using principal-component analysis based on the top 2,000 highly variable genes by the ‘FindVariableGenes’ function in Seurat. Cell clusters on UMAP were identified using the ‘FindClusters’ function of Seurat.

Pseudotime analysis

We conducted pseudotime analysis using Monocle3 (ref. 67). Highly variable genes identified by ‘FindVariableGenes’ function of Seurat were used to estimate similarities between cells. Pseudotimes were calculated according to similarities in the transcriptome or epigenome to reflect cell relationships during differentiation. Further, cells were ordered by pseudotime as shown in Figs. 2d and 3d.

Genomic annotation analysis

Genes near to regulatory elements were selected with the annotatePeaks.pl function in Homer68. GO analysis for biological processes was performed using Metascape (http://metascape.org)69, with Benjamini–Hochberg P value correction70.

Data downsampling

The sampling rates were set at predetermined values, including 80% to 20%. Each sampling rate corresponded to an independent downsampling process. Utilizing the SAMtools view command (with -s specifying the sampling rate and -b for generating BAM format files), a proportionate subset of reads was randomly selected from the original BAM files to create new BAM files for each sampling rate.

Cluster error rate calculation

A cisTopic object was created from the downsampled BAM files and genomic region files in BED format using the createcisTopicObjectFromBAM function. Collapsed Gibbs Sampling models were then executed to explore the optimal model selection among cistopics. Model selection was assessed based on the log-likelihood iteration plot. For the finalized model, dimensionality reduction was conducted using the UMAP method, which subsequently formed the basis for clustering analysis of cells. Hierarchical clustering was performed using the fastcluster package, categorizing cells into two groups based on their UMAP coordinates. This error rate was defined as the proportion of cells incorrectly grouped relative to the total number of cells. Specifically, clustering outcomes were compared against known cell types, with any incongruent allocations considered erroneous.

Statistics and reproducibility

All EpiChem experiments were independently performed at least twice. The Seurat package implemented in R was used to identify variable genes by calculating the variance and finding the differentially expressed. The Monocle 3 package was used to calculate pseudotime as described above. Other statistical analyses for EpiChem were mainly performed between groups of peaks or segments of the genome. The statistical tests are indicated in the figure legends.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.