{"id":593620,"date":"2023-01-03T06:49:28","date_gmt":"2023-01-03T12:49:28","guid":{"rendered":"https:\/\/news.sellorbuyhomefast.com\/index.php\/2023\/01\/03\/scchix-seq-infers-dynamic-relationships-between-histone-modifications-in-single-cells\/"},"modified":"2023-01-03T06:49:28","modified_gmt":"2023-01-03T12:49:28","slug":"scchix-seq-infers-dynamic-relationships-between-histone-modifications-in-single-cells","status":"publish","type":"post","link":"https:\/\/newsycanuse.com\/index.php\/2023\/01\/03\/scchix-seq-infers-dynamic-relationships-between-histone-modifications-in-single-cells\/","title":{"rendered":"scChIX-seq infers dynamic relationships between histone modifications in single cells"},"content":{"rendered":"
\n
\n

Main<\/h2>\n
\n

Gene expression in animals relies on epigenetic marks such as histone modifications to regulate the accessibility and function of the genome in different cell types1<\/a><\/sup>. Large-scale efforts characterizing different histone modifications in a variety of cell populations commonly use chromatin immunoprecipitation followed by sequencing (ChIP\u2013seq)2<\/a>,3<\/a>,4<\/a>,5<\/a>,6<\/a>,7<\/a>,8<\/a><\/sup>. Alternative strategies to ChIP\u2013seq based on enzyme tethering (chromatin immunocleavage, ChIC) have reduced the background signal in profiling the epigenome9<\/a><\/sup>, and have enabled single-cell profiling of histone modifications8<\/a>,10<\/a>,11<\/a>,12<\/a>,13<\/a>,14<\/a>,15<\/a>,16<\/a>,17<\/a>,18<\/a>,19<\/a><\/sup>. Tethering strategies involve incubating cells with an antibody against a histone modification of interest, which then tethers either protein A-MNase10<\/a>,12<\/a>,18<\/a>,19<\/a><\/sup> or protein A-Tn511<\/a>,13<\/a>,14<\/a>,15<\/a>,16<\/a>,17<\/a><\/sup> fusion protein to generate targeted DNA fragments in single cells. However, most experimental techniques to map single-cell histone modifications are limited to only one histone modification per single cell.<\/p>\n

We present an integrated experimental and computational framework for multiplexing histone modifications in single cells. To profile two histone modifications in single cells (Fig. 1a<\/a>), we first generate three genome-wide sortChIC18<\/a><\/sup> datasets: two datasets by incubating cells with one of the two histone modification antibodies separately (single-incubated; Fig. 1b<\/a>), and the third by incubating cells with both histone modification antibodies together (double-incubated; Fig. 1b<\/a>). We then use our two single-incubated datasets as training data to generate the possible pairs of genome-wide histone modification profiles that, when added together, fit to a single-cell profile from the double-incubated dataset (Fig. 1c<\/a>). For each double-incubated cell, we then deconvolve the multiplexed data by probabilistically assigning each fragment back to their respective histone modification.<\/p>\n

\n
Fig. 1: Overview of the scChIX-seq method.<\/b><\/figcaption>
\n
\"figure<\/picture><\/a><\/div>\n

a<\/b>, Chromatin regulation of different cell types (different colored cells) is regulated in part through several histone modifications (two histone modifications shown as an example). b<\/b>, scChIX-seq uses three sortChIC antibody incubation conditions: two conditions each target a single histone modification (single-incubated) only and the third condition targets both histone modifications simultaneously (double-incubated). c<\/b>, Schematic of scChIX-seq for deconvolving multiplexed histone modifications. The two single-incubated sortChIC datasets (one targeting an orange histone modification, the other a blue modification, each modification reveals three clusters) are training data to define the possible pairs of histone modification distributions that can be combined to generate a hypothetical double-incubated cell. For each observed double-incubated cell, we then assign the cell to the most probable pair of cell states, one from each histone modification. We then probabilistically assign each pA-MNase cut into their respective histone modification. Cartoons represent genome-wide distribution of histone modification signals in different modifications and cell types; x<\/i> axes represent genomic distance, and vertical ticks are arbitrary distance markers. d<\/b>, Label transfer allows joint analysis of two single-incubated sortChIC datasets targeting functionally distinct histone modifications. Information derived from one histone modification, such as cell types, histone mark levels and pseudotime, can be transferred to another histone modification using the double-incubated cells as a link. e<\/b>, Simulation study shows that scChIX-seq can unbiasedly assign reads to each mark regardless of the amount of overlap there is between the two marks across the genome. x<\/i> axis of cartoon genome-wide distributions (middle-left) is genomic distance. Right: ground truth probabilities versus inferred probabilities from scChIX. p<\/i> is the expected fraction of double-incubated reads in a genomic locus that belongs to mark 1. (hat{p})<\/span> is the estimate of the probability; n<\/i>\u2009=\u2009101 simulation datapoints spread evenly between 0 and 1 inclusive. Error bars are 95% CI, centers are the mean.<\/p>\n<\/div>\n

Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n

scChIX-seq links single-cell maps of different histone modifications, revealing relationships between histone modifications in single cells. In these linked maps, information derived from one chromatin state, such as cell types, histone mark levels and pseudotimes, can transfer to another chromatin state (Fig. 1d<\/a>), unlocking joint analysis of several histone modifications in single cells. We first validated scChIX-seq using simulation, purified blood cell types and whole bone marrow. We then applied scChIX-seq to two complex biological systems, one in mouse organogenesis to uncover orthogonal dynamics in H3K36me3 and H3K9me3, and the other in macrophage in vitro differentiation to reveal coordinated dynamics between H3K4me1 and H3K36me3.<\/p>\n<\/div>\n<\/div>\n

\n

Results<\/h2>\n
\n

Benchmarking across histone modification relationships<\/h3>\n

To test whether scChIX-seq is accurate for histone modification patterns that are mutually exclusive as well as highly overlapping, we apply scChIX-seq to simulated single-cell data with known amounts of overlap to benchmark our method across different overlapping patterns between histone modifications. We simulate single-cell histone modification data by modifying simATAC20<\/a><\/sup> to generate sparse count data from different overlapping patterns from the same cell (Fig. 1e<\/a> and Extended Data Fig. 1a,b<\/a>; Methods<\/a>). Our simulations span three scenarios to cover varying degrees of overlapping patterns (Extended Data Fig. 1c<\/a>). (1) Mutually exclusive scenario with only 1% of loci overlapping. (2) Intermediate scenario with 50% of loci overlapping. (3) Correlated scenario with 99% of loci overlapping. In these simulations, we provide a ground truth parameter p<\/i> for each genomic locus and then estimate this parameter using our statistical framework to assess the uncertainty in our inferences. Here, p<\/i> is the expected fraction of double-incubated reads in a locus that belongs to a reference histone modification (that is, p<\/i>\u2009=\u20090.5 if locus is exactly overlapping, p<\/i>\u2009=\u20091 or 0 if locus is exactly mutually exclusive). Applying scChIX-seq to each scenario, we find that the distribution of our estimates (hat{p})<\/span> across all loci are comparable with the ground truth distribution of p<\/i> (Extended Data Fig. 1c,d<\/a>). Furthermore, scChIX-seq accurately recovers the different cell types underlying the simulated data, and links the two histone modification landscapes into a joint uniform manifold approximation and projection (UMAP) (Extended Data Fig. 1e<\/a>). Summarizing the three scenarios, scChIX-seq can estimate p<\/i> accurately for all degrees of overlap, with confidence intervals (CI) better than (hat{p}pm 0.05)<\/span> (Fig. 1e<\/a> (right) and Extended Data Fig. 1f<\/a>). Our simulation study confirms that scChIX-seq is accurate in inferring several histone modifications in single cells in both mutually exclusive as well as overlapping histone modification patterns.<\/p>\n

Validating with ground truth data from purified cell types<\/h3>\n

To validate our method experimentally, we generate a ground truth sortChIC dataset by purifying three known cell types from mouse bone marrow: B cells, granulocytes and natural killer (NK) cells, using fluorescence-activated cell sorting (FACS) and applying scChIX-seq (Methods<\/a>). Of note, the sortChIC method is designed to integrate FACS with histone modification mapping18<\/a><\/sup>, so we can enrich for a cell type and map histone modifications in one workflow. We split bone marrow cells into three technical batches: one batch incubated with anti-H3K27me3 antibody alone (single-incubated), one with anti-H3K9me3 alone (single-incubated) and the third with both anti-H3K27me3 and anti-H3K9me3 antibodies together (double-incubated, H3K27me3+H3K9me3). We then sorted cells into 384-well plates, each plate containing all three cell types, and generate targeted cut fragments (Extended Data Fig. 2a,b<\/a>). We chose H3K27me3 and H3K9me3 because they have been shown to have a mutually exclusive relationship21<\/a><\/sup>, allowing us to verify whether we can infer the correct cell type as well as the generally mutually exclusive relationship. Of note, although H3K27me3 and H3K9me3 are known to be nonoverlapping, it is unclear how this relationship precisely changes to make cell-type-specific patterns at different loci, and therefore modeling the two relationships is still needed to accurately infer the two chromatin profiles in individual cells.<\/p>\n

From the double-incubated data alone, we would not know which cut fragments correspond to H3K27me3 and which to H3K9me3, but would observe only a superposition of the two profiles. We therefore used the single-incubated sortChIC data to train a statistical model of how cells from the same cell type combine their H3K27me3 and H3K9me3 profiles to generate double-incubated cut fragments. This model was then used to deconvolve the single-cell multiplexed signal into their respective histone modifications (Methods<\/a>).<\/p>\n

To learn an interpretable latent space for H3K27me3 and H3K9me3, we applied latent Dirichlet allocation (LDA)22<\/a>,23<\/a><\/sup> to the single-incubated H3K27me3 and H3K9me3 datasets, which factorizes count matrices based on a multinomial model (Methods<\/a>). (Extended Data Fig. 2c,d<\/a>). LDA learns cell-type-specific vectors of probabilities. These parameters model the probability that a cut fragment would fall into a specific genomic region. These probabilities can therefore be interpreted as genome-wide histone modification distributions that depend on cell type, and each cell generates a high-dimensional sparse count vector with n<\/i> total fragments by drawing n<\/i> independent trials from these multinomial distributions.<\/p>\n

Demultiplexing the double-incubated data involves two steps. First, we used the training data to infer which genome-wide H3K27me3 distribution was added to which H3K9me3 distribution to generate a linear combination of two distributions (H3K27me3+H3K9me3). Second, we probabilistically assigned each double-incubated cut fragment to either H3K27me3 or H3K9me3, given that we know the underlying linear combination of the two profiles.<\/p>\n

The deconvolved H3K27me3+H3K9me3 data generated two sets of cuts for each cell: one set coming from H3K27me3 and the other from H3K9me3. We projected the two sets of cuts onto the H3K27me3 or H3K9me3 latent space (learned from LDA), respectively (Fig. 2a<\/a>). Since each deconvolved cell has a set of cuts in H3K27me3 and H3K9me3 simultaneously, we can link the UMAPs together, creating a joint chromatin regulation space (Fig. 2a<\/a>).<\/p>\n

\n
Fig. 2: scChIX-seq accurately deconvolves multiplexed histone modifications in single cells.<\/b><\/figcaption>
\n
\"figure<\/picture><\/a><\/div>\n

a<\/b>, UMAP representation of the H3K27me3 (n<\/i>\u2009=\u2009367) and H3K9me3 (n<\/i>\u2009=\u2009376) histone modification space derived from the two single-incubated datasets (right two panels), and the H3K27me3+H3K9me3 space (left panel, n<\/i>\u2009=\u2009290) derived from the double-incubated data. Cells are colored by their ground truth cell-type labels. The cells in the H3K27me3- and H3K9me3-only space have unmixed double-incubated cells whose deconvolved signal has been projected onto their respective UMAPs. Lines connecting across datasets connect where each double-incubated cell is located in each of the three histone modification space. b<\/b>, Matrix summarizing the cluster pair that scChIX-seq selected for each double-incubated cell. Cells along the diagonal are predicted to be B cells, granulocytes and NK cells, respectively. Cells in the off-diagonal are false negatives. Barplots summarizing FDR, sensitivity and specificity of assigning each cell type (right). c<\/b>, Zoom-in coverage plot and single-cell cut fragments in B cells of mixed (H3K27me3+H3K9me3, gray bars), unmixed (H3K27me3 and H3K9me3, orange and blue bars). Positions of cut fragments are shown for four single cells (single cells A, B, C and D) for H3K27me3+H3K9me3 signal (gray ticks) as well as their unmixed outputs (orange and blue ticks). Circled reads and arrow highlight examples of cut fragments being assigned to either H3K27me3 (orange) or H3K9me3 (blue). d<\/b>, Zoom-out of the Serpinb5<\/i> locus. Cut fragments from H3K27me3+H3K9me3 are colored based on whether they have been assigned to H3K27me3 (orange) or H3K9me3 (blue). Ground truth coverage are single-incubated sortChIC data targeting H3K27me3 (orange) and H3K9me3 (blue). e<\/b>, Heatmap of probabilities p<\/i> of assigning reads to H3K27me3 (p<\/i>\u2009=\u20091, red) or H3K9me3 (p<\/i>\u2009=\u20090, blue) around the Bcl2<\/i> locus. Rows are single cells (ordered by predicted cell type), columns are genomic regions (50\u2009kb bins). Transitions between H3K9me3- and H3K27me3-marked chromatin states are independent of cell type. f<\/b>, Same as e<\/b> but at the Crim1<\/i> locus, where transitions from H3K9me3 to H3K27me3 (blue to red) are cell-type specific.<\/p>\n<\/div>\n

Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n

The double- and single-incubated cells in the H3K27me3 and H3K9me3 UMAPs intermingle, suggesting that the model accurately assigns cut fragments to their respective histone modification (Extended Data Fig. 2e,f<\/a>). Comparing the H3K27me3 deconvolved pseudobulk signal with our ground truth single-incubated pseudobulk shows high correlation for the expected cell type, and lower for the other two cell types (Extended Data Fig. 2g<\/a>). The H3K9me3 deconvolved pseudobulk signal also shows highest correlation with the expected cell type, with lower correlation from other cell types (Extended Data Fig. 2h<\/a>). Finally, we compared the fragments per cell obtained from scChIX-seq versus multi-CUT&TAG24<\/a><\/sup>, and found that scChIX-seq achieves higher sensitivity than multi-CUT&TAG (Extended Data Fig. 2i). Overall, our ground truth dataset demonstrates that scChIX-seq is accurate and sensitive in assigning cut fragments to their respective histone modification.<\/p>\n

To quantify the accuracy of scChIX-seq in selecting the correct H3K27me3-H3K9me3 cluster pair to mix together, we color each cell by its ground truth label and plot its inferred H3K27me3-H3K9me3 pair on a two-dimensional (2D) grid (Fig. 2b<\/a>, left). The false discovery rates (FDRs) of scChIX-seq predicting B cells, granulocytes or NK cells are 10%, 3% and 1%, respectively (Fig. 2b<\/a>, right). Similarly, scChIX-seq has high specificity and sensitivity in inferring the correct cluster pairs (Fig. 2b<\/a>, right).<\/p>\n

Next, scChIX-seq assigns each double-incubated cut fragment to either H3K27me3 or H3K9me3 (Fig. 2c<\/a>; Methods<\/a>). The deconvolved B cell repressive landscapes correspond with their respective ground truth, exemplified in the Bcl2<\/i> (Fig. 2d<\/a>) and Crim1<\/i> (Extended Data Fig. 3a<\/a>) locus. We also find cell-type-specific signal in H3K27me3 (Extended Data Fig. 3b<\/a>) and H3K9me3 signal (Extended Data Fig. 3c<\/a>).<\/p>\n

Our model infers p<\/i>, the expected fraction of double-incubated fragments at a locus that belongs to H3K27me3. That is, p<\/i>\u2009=\u20090 if all fragments belong to H3K9me3 and p<\/i>\u2009=\u20091 if they all belong to H3K27me3. Plotting these probabilities across all loci reveals a bimodal distribution with peaks near 0 and 1 (Extended Data Fig. 3d<\/a>). Classifying these loci as H3K9me3-specific (P<\/i>\u2009<\u20090.5) or H3K27me3-specific (P<\/i>\u2009\u2265\u20090.5), we compare the GC content and distance to transcription start site (TSS) of the two classes of loci (Extended Data Fig. 3e,f<\/a>). We find H3K9me3-specific regions to have lower GC content and increased distance from TSSs compared with H3K27me3-specific regions. Of note, we observe this difference across all three cell types, suggesting that GC-poor and gene-poor regions of the genome is a general feature of H3K9me3-specific regions21<\/a><\/sup>.<\/p>\n

Summarizing these probabilities in single cells along the genome as a heatmap, the Bcl2<\/i> locus reveals the mutual exclusive relationship between H3K27me3 and H3K9me3, where the chromatin state is predominantly H3K9me3, then switches to H3K27me3, and then switches back to H3K9me3 (Fig. 2e<\/a>). For Bcl2<\/i>, these transitions occur at the same location independent of the cell type. However, we also find that these transitions can be cell-type specific, as exemplified by the Crim1<\/i> locus (Fig. 2f<\/a>), where the H3K27me3 region extends further upstream of Crim1<\/i> in NK cells compared with B cells and granulocytes. Our ground truth experiment demonstrates that scChIX-seq can accurately map two histone modifications in single cells, and the inferred probabilities can be biologically interpreted as relationships between the two histone modifications in single cells.<\/p>\n

scChIX-seq reveals H3K4me1\/H3K27me3 relationships in bone marrow<\/h3>\n

We next apply scChIX-seq to integrate active (H3K4me1) and repressive (H3K27me3) chromatin states in a complex mixture of cells by sampling mouse bone marrow (Extended Data Fig. 4a,b<\/a>). We use scChIX-seq to transfer labels and link UMAPs between active and repressive histone modifications (Fig. 3a,b<\/a>) to perform a joint analysis of the two marks.<\/p>\n

\n
Fig. 3: scChIX-seq enables joint analysis of distinct histone modifications in single cells.<\/b><\/figcaption>
\n
\"figure<\/picture><\/a><\/div>\n

a<\/b>, UMAP of sortChIC signal of H3K4me1 in bone marrow (n<\/i>\u2009=\u2009639 cells). Clusters are colored by cell type. Latent space calculated using LDA with 50\u2009kb bins. b<\/b>, UMAP of sortChIC signal of H3K27me3 in whole bone marrow (n<\/i>\u2009=\u2009517 cells). Cell types in H3K27me3 are inferred by transferring labels from H3K4me1. c<\/b>, H3K4me1 and H3K27me3 UMAPs linked together by deconvolved double-incubated cells (n<\/i>\u2009=\u20091,711 cells). H3K4me1 and H3K27me3 portions of the double-incubated cells are projected onto their respective UMAPs. Lines connect where the active signal and the corresponding repressive signal are located for each double-incubated cell. DC, dendritic cells; pDC, plasmacytoid dendritic cells. d<\/b>, Heatmap showing probability of assigning a read in a region to either H3K27me3 or H3K4me1 at 5\u2009kb resolution. Heatmap shows the Igk<\/i> locus for pro-B versus B cells. Rows are single cells, columns are 5\u2009kb genomic regions. Blue represents regions where cut fragments are probably coming from H3K27me3, while red represents regions where cut fragments are probably coming from H3K4me1.<\/p>\n<\/div>\n

Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n

To define cell types from the H3K4me1 sortChIC data, we ranked the top 150 genes associated with different clusters from sortChIC and used a publicly available scRNA-seq dataset to compare mRNA abundances of cluster-specific genes across different blood cell types25<\/a><\/sup> (Extended Data Fig. 4c<\/a>). scChIX-seq takes each H3K4me1+H3K27me3 cell and infers the most probable cluster pair (one from H3K4me1, the other from H3K27me3), which systematically transfers cell-type labels defined from H3K4me1 onto the H3K27me3 data (Extended Data Fig. 4d<\/a>). We find that a small minority of double-incubated cells have low-confidence cluster pair predictions. Plotting the cluster pairs onto the H3K4me1+H3K27me3 UMAP confirms that the single-cell assignment produces precise clusters where neighboring cells are probably assigned to the same pair. Low-confidence predictions arise from cells that border between clusters (Extended Data Fig. 4e<\/a>), which we remove from further analysis. Overall, scChIX-seq allows systematic transfer of cell-type labels from one histone modification to another.<\/p>\n

We next deconvolve the double-incubated cells into their respective histone modification. The UMAPs from H3K4me1 and H3K27me3 show that single-incubated and deconvolved single cells intermingle, suggesting that deconvolution does not produce batch effects (Extended Data Fig. 4f,g<\/a>). The deconvolved single cells provide anchors to systematically link one histone modification with another (Fig. 3c<\/a>). To validate the predicted cell types in both the single and deconvolved datasets, we compared with data from cell types purified by FACS. For H3K4me1 clusters, we compared with publicly available ChIP\u2013seq5<\/a><\/sup>. Pearson correlation between ChIP\u2013seq of B cells, erythroids, granulocytes and NK cells versus sortChIC from single- and double-incubated cells is highest for the predicted cell type (Extended Data Fig. 5a\u2013d<\/a>). Although single-incubated cells have higher correlation with ChIP\u2013seq reference data than deconvolved cells for the matched cell type, the deconvolved cells of the matched cell type consistently had higher correlation with ChIP\u2013seq than unmatched cell types. For H3K27me3 clusters, we used our ground truth sortChIC data purified from FACS. Pearson correlation of sortChIC signal between FACS-sorted B cells, granulocytes and NK cells versus pseudobulks derived from whole bone marrow is highest for the predicted cell type (Extended Data Fig. 5e\u2013g<\/a>).<\/p>\n

Classifying these loci as H3K27me3-specific or H3K4me1-specific using a cluster-specific cutoff for p<\/i> (Extended Data Fig. 5h<\/a>), we again compare the GC content and distance to TSS of the two classes of loci. We find that H3K4me1-marked regions tend to be closer to TSSs compared with H3K27me3 (Extended Data Fig. 5i<\/a>), and that GC content is higher in H3K27me3-specific compared with H3K4me1-specific regions (Extended Data Fig. 5j<\/a>). The increase in GC content for H3K27me3-marked regions is consistent with previous studies showing that GC-rich elements in transcriptionally inactive regions can recruit PRC2 (ref. 26<\/a><\/sup>).<\/p>\n

We use the joint landscape to reveal active and repressive histone modification dynamics within cell types. To find differences in chromatin regulation between pro-B cells versus B cells, we select only pro-B or B cells and recluster the cells in both H3K4me1 and H3K27me3 separately (Extended Data Fig. 6a,b<\/a>). With multimodal data, we can transfer cell-type-specific H3K4me1 signal onto the H3K27me3 UMAP to distinguish pro-B and B cells with more confidence. Using pro-B cell-specific genes, Pax5<\/i> (ref. 27<\/a><\/sup>) and Pten<\/i>28<\/a><\/sup>, we project the H3K4me1 signal at loci overlapping these genes onto both H3K4me1 and H3K27me3 landscapes, confirming a subset of pro-B cells within the B cell population (Extended Data Fig. 6c<\/a>). Similarly, we use marker genes associated with more differentiated B cells, such as Irf4<\/i> (ref. 27<\/a><\/sup>), Igkv3-2<\/i> locus29<\/a><\/sup> and Cd72<\/i> (ref. 30<\/a><\/sup>) to confirm a more differentiated B cell population (Extended Data Fig. 6d<\/a>). Plotting the heatmap of H3K4me1-H3K27me3 assignment probabilities at the IgK<\/i> locus reveals that the chromatin state is repressed in pro-B cells but becomes activated in B cells (Fig. 3d<\/a>), consistent with the progressive activation of the chromatin state during B cell development29<\/a><\/sup>.<\/p>\n

Next, we recluster neutrophils to analyze differences in chromatin regulation along pseudotime (Extended Data Fig. 7a<\/a>). Reclustering neutrophils in H3K27me3 reveals a shared pseudotime trajectory that varies smoothly between neutrophils in both the H3K27me3 and H3K4me1 landscapes. H3K4me1 levels at the Retnlg<\/i> locus\u2014a marker gene for mature neutrophils31<\/a><\/sup>\u2014increases along pseudotime, while H3K27me3 levels decreases (Extended Data Fig. 7b<\/a>). The H3K27me3 gene loadings associated with pseudotime consists of a module of Hox<\/i> and other developmental genes (Extended Data Fig. 7c\u2013e<\/a>). Of note, these genes have low levels of mRNA abundances in neutrophils (Extended Data Fig. 7f<\/a>), suggesting that this module is transcriptionally silent. At a locus overlapping the Hoxa<\/i> locus, we find that H3K27me3 was highly marked while H3K4me1 was lowly marked across all neutrophils. Along pseudotime, H3K27me3 increases further, while H3K4me1 decreases further (Extended Data Fig. 7c<\/a>). Our pseudotime analysis suggests that dynamics in histone modifications can occur even in regions associated with low-expressed genes.<\/p>\n

H3K36me3\/H3K9me3 relationships during mouse organogenesis<\/h3>\n

To demonstrate the method in more complex biological scenarios, we applied scChIX-seq during mouse organogenesis (E9.5 to E11.5) to study H3K36me3 and H3K9me3 dynamics at single-cell resolution (Fig. 4a<\/a> and Extended Data Fig. 8a,b<\/a>). We took the top 250 cluster-specific bins from the H3K36me3 data to identify cell types (Methods<\/a>). These loci associate with gene bodies of cell-type-specific genes. For example, we find H3K36me3 signal around genes enriched in specific cell types, such as erythroids (Sptb<\/i>)32<\/a><\/sup>, white blood cells (Lcp2<\/i> (ref. 33<\/a><\/sup>), endothelial cells (Emcn<\/i>)34<\/a><\/sup>, neural tube (Rfx4<\/i>)35<\/a><\/sup>, neurons (Elavl4<\/i>)36<\/a><\/sup>, Schwann precursors (Cdh6<\/i>)37<\/a><\/sup>, epithelial cells (Grhl2<\/i>)38<\/a><\/sup>, mesenchymal progenitors (Prx1<\/i>)39<\/a><\/sup> and cardiomyocytes (Gata6<\/i>, Tpm1<\/i>)40<\/a>,41<\/a><\/sup> (Extended Data Fig. 8c\u2013l<\/a>).<\/p>\n

\n
Fig. 4: Applying scChIX-seq to mouse organogenesis reveals shared heterchromatin landscapes and cell-type-specific differences in H3K36me3:H3K9me3 ratios.<\/b><\/figcaption>
\n
\"figure<\/picture><\/a><\/div>\n

a<\/b>, Schematic of mouse organogenesis experiment to study H3K36me3 and H3K9me3 in single cells. b<\/b>, Joint UMAP of mouse organogenesis after deconvolution from scChIX-seq (n<\/i>\u2009=\u20092,911 H3K36me3 cells, n<\/i>\u2009=\u20092,166 H3K9me3 cells). c<\/b>, Assignment of several H3K36me3 cell types to one H3K9me3 cluster. The H3K36me3 (columns) and H3K9me3 (rows) label for each double-incubated cells (n<\/i>\u2009=\u20091,197 cells) are plotted onto a matrix to H3K36me3 cell types to H3K9me3 clusters. Cells are colored by their cell-type label as in b<\/b>. d<\/b>, Subclustering of nonblood cells for H3K9me3, represented as a UMAP. Arrow denotes a pseudotime axis. Pseudotime defined as the first PC of the 2D UMAP. e<\/b>, Joint UMAP of deconvolved double-incubated cells (n<\/i>\u2009=\u20091,197 cells), colored by the log ratio of number of H3K36me3 cuts versus number of H3K9me3 cuts. f<\/b>, Boxplot of H3K36me3:H3K9me3 ratio across cell types. Number of double-incubated cells for each cell type: n<\/i>\u2009=\u2009163 erythroid, n<\/i>\u2009=\u200917 white blood cells, n<\/i>\u2009=\u200924 endothelial, n<\/i>\u2009=\u2009136 neural tube progenitors, n<\/i>\u2009=\u2009197 neurons, n<\/i>\u2009=\u200946 Schwann cell precursors, n<\/i>\u2009=\u200973 epithelial, n<\/i>\u2009=\u2009458 mesenchymal progenitors and n<\/i>\u2009=\u200983 cardiomyocytes. Boxplots show 25th percentile, median and 75th percentile, with the whiskers spanning 97% of the data.<\/p>\n<\/div>\n

Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n

To uncover whether distinct H3K36me3 cell types could share common H3K9me3 landscapes, we deconvolved the H3K36me3 + H3K9me3 cells and projected each cell to both landscapes (Fig. 4b<\/a>). scChIX-seq reveals that erythroid and white blood cells have both distinct active chromatin and heterochromatin, but the other nonblood cell types show similar heterochromatin distribution. Assigning each double-incubated cell to a H3K36me3 and H3K9me3 cluster confirms that cells with distinct H3K36me3 can share the same H3K9me3 cluster (Fig. 4c<\/a>). Of note, the variable genes that show cell-type-specific differences in both active chromatin and publicly available mRNA abundances42<\/a><\/sup> (Extended Data Fig. 9a,b<\/a>) have low signal across cell types in H3K9me3 (Extended Data Fig. 9c<\/a>), suggesting that using conventional marker genes from RNA-seq would not reveal cell-type differences in H3K9me3.<\/p>\n

Differential expression across the three H3K9me3 clusters reveals cluster-specific repressed loci (Extended Data Fig. 9d<\/a>), with the largest effect coming from erythroid-specific regions. These erythroid-repressed regions are associated with decreased mRNA abundances (Extended Data Fig. 9e\u2013g<\/a>). Subsetting the data and running LDA on only nonblood cells in H3K9me3, we find that H3K9me3 varies over organogenesis stages (Fig. 4d<\/a>), suggesting that heterochromatin differences are stronger across organogenesis stages than between cell types.<\/p>\n

Because the double-incubated cells have cut fragments associated with both histone modifications, we hypothesized that the deconvolved data could precisely quantify the ratio between the two histone modifications, and how this ratio changes across cell types. Counting total reads from single-incubated data would lead to large cell-to-cell technical variability because counts per cell can span several orders of magnitude. However, comparing the counts of the two histone modification in the same cell could overcome this technical variability. We therefore asked whether the global ratio of H3K36me3 and H3K9me3 in individual cells varies. Plotting the ratio of H3K36me3 and H3K9me3 reveals that most cells have comparable ratios, but that erythroid cells have lower ratios than other cell types (Fig. 4e,f<\/a>). This lower ratio is consistent with mass spectrometry studies showing a global decrease in H3K36me3 but no change in H3K9me3 during erythroid maturation43<\/a><\/sup>. Of note, inferring this global change without scChIX-seq, such as by counting total unique fragments from single-incubation data, is challenging due to the large variability in total counts across cells and the inability to distinguish cell types in certain H3K9me3 clusters (Extended Data Fig. 9h,i<\/a>).<\/p>\n

In sum, applying scChIX-seq to H3K36me3 and H3K9me3 during organogenesis reveals unique insights from multimodal analysis. The complex relationships between the two histone modifications as well as their global changes would not have been elucidated by analyzing single-incubated data alone.<\/p>\n

Mark-specific pseudotimes and chromatin velocity<\/h3>\n

Finally, we applied scChIX-seq to study the dynamic relationships between two active histone modifications, H3K4me1 and H3K36me3, over an in vitro differentiation timecourse. We sorted blood progenitors from mouse bone marrow, added macrophage colony-stimulating factor (MCSF) and collected cells over 7\u2009days (Fig. 5a<\/a> and Extended Data Fig. 10a,b<\/a>; Methods<\/a>). We incubated cells with either H3K4me1, H3K36me3 or both H3K4me1 and H3K36me3, then performed scChIX-seq.<\/p>\n

\n
Fig. 5: Applying scChIX-seq to two active marks reveals chromatin velocity during in vitro macrophage differentiation.<\/b><\/figcaption>
\n
\"figure<\/picture><\/a><\/div>\n

a<\/b>, Schematic of mouse macrophage in vitro differentiation timecourse experiment to study H3K4me1 and H3K36me3 in single cells. b<\/b>, Heatmap of histone modification signal on the bodies of dynamic genes over pseudotime. Rows are gene bodies and columns are single-incubated cells ordered along pseudotime. Color labels of columns are days from which the cells were recovered during the timecourse. c<\/b>, Boxplots of pseudotime estimates of single-incubated cells along the timecourse. Number of cells per day for H3K4me1: n<\/i>\u2009=\u200958 day 0, n<\/i>\u2009=\u2009148 day 1, n<\/i>\u2009=\u2009249 day 2, n<\/i>\u2009=\u2009350 day 3, n<\/i>\u2009=\u2009369 day 4, n<\/i>\u2009=\u2009383 day 5, n<\/i>\u2009=\u2009488 day 6, n<\/i>\u2009=\u2009519 day 7. For H3K36me3: n<\/i>\u2009=\u200942 day 0, n<\/i>\u2009=\u2009125 day 1, n<\/i>\u2009=\u2009176 day 2, n<\/i>\u2009=\u2009301 day 3, n<\/i>\u2009=\u2009384 day 4, n<\/i>\u2009=\u2009366 day 5, n<\/i>\u2009=\u2009522 day 6, n<\/i>\u2009=\u2009567 day 7. Boxplots show 25th percentile, median and 75th percentile, with the whiskers spanning 97% of the data. d<\/b>, Estimate of the average difference of pseudotime from one day to the next. Error bars indicate 95% CI, calculated by a linear model of the pseudotime differences between days. Statistics derived from number of cells indicated in c<\/b>. e<\/b>, Estimates of two different pseudotimes from a single cell. Error bars are 95% CI of the estimates. Each point is a double-incubated cell. f<\/b>, Joint UMAP of H3K4me1 and H3K36me3 from scChIX-seq, lines connect single cells with multimodal information. g<\/b>, Chromatin velocity estimates of an upregulated gene (above) and a downregulated gene (below). Red curve is the exponential relaxation fit according to the solution of the first-order differentiation equation. h<\/b>, High-dimensional chromatin velocities of dynamic genes projected onto PCs 1 and 2. Vector field estimated by smoothing across nearest neighbors of cells (Methods<\/a>).<\/p>\n<\/div>\n

Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n

Genome tracks of H3K4me1 and H3K36me3 signal for each day shows upregulation of macrophage-specific genes, such as Mertk<\/i>44<\/a><\/sup> (Extended Data Fig. 10c<\/a>). Heatmap of H3K4me1 and H3K36me3 dynamics at gene bodies along pseudotime reveals that the two histone modifications up- and downregulate genes with different dynamics. H3K36me3 shows a gradual up- or downregulation of signal while H3K4me1 reaches a new steady state earlier along pseudotime (Fig. 5b<\/a>). Summarizing log2<\/sub> fold change of the two histone modifications genome-wide, we find that dynamics in H3K36me3 are often larger than in H3K4me1 (Extended Data Fig. 10d<\/a>). Comparing pseudotime progression with day of sample collection shows that changes in H3K4me1 peak at day\u20092 and then increases progressively over the day while H3K36me3 dynamics peak around day\u20093 and 4 before relaxing towards steady state (Fig. 5c<\/a>). The time of the largest change in H3K4me1 dynamics occurs 1\u2009day before H3K36me3 (Fig. 5d<\/a>), suggesting that global changes in H3K4me1 precede changes in H3K36me3. Summarizing at the genome-wide level, UMAPs of H3K4me1 and H3K36me3 of single-incubated cells show that both active marks move progressively towards a macrophage state during the timecourse (Fig. 5e<\/a>).<\/p>\n

Using continuous pseudotime of H3K4me1 and H3K36me3 as our training data (Methods<\/a>), for both H3K4me1 and H3K36me3 we infer where along pseudotime each double-incubated cell came from. Plotting the inferred pseudotimes of each mark for each cell uncovers the dynamic relationships between the two marks (Fig. 5e<\/a>). H3K4me1 pseudotime initially progresses while H3K36me3 remains relatively unchanged. As H3K4me1 pseudotime approaches 0.5, H3K36me3 then progresses rapidly towards 1. This sigmoidal-like relationship between H3K4me1 versus H3K36me3 pseudotime progression is consistent with H3K4me1 dynamics occurring before H3K36me3. Finally, we used this inferred pseudotime information to project the deconvolved cells onto the H3K4me1 and H3K36me3 UMAPs. Both UMAPs showed that the single-incubated and deconvolved cells intermingle with each other, suggesting that deconvolution was successful (Extended Data Fig. 10e,f<\/a>). Using the deconvolved cells as anchors, we then linked the two histone modification maps together (Fig. 5f<\/a>).<\/p>\n

Since we observed that H3K4me1 dynamics occur before H3K36me3, we asked whether we could model the H3K36me3 dynamics as a first-order differential equation analogous to RNA velocity45<\/a><\/sup> (Fig. 5g<\/a>, top; Methods<\/a>). Since our data come from a timecourse, we directly fitted the exponential curves for dynamic genes along pseudotime for H3K36me3 (Extended Data Fig. 10g<\/a>), which avoids making steady-state assumptions and leverages information from both single-incubated and deconvolved cells across histone modifications. The distribution of inferred rate constants from the exponential fit show a median of approximately 2.3 per pseudotime (Extended Data Fig. 10h<\/a>). These rate constants were then used to predict the H3K36me3 levels for each cell over small pseudotime steps (\u0394t<\/i> = 0.02; Fig. 5g<\/a>). Finally, summarizing the predictions of dynamic genes, we projected the high-dimensional velocity vectors onto the first two principal components (PCs). From the chromatin velocity summary, we found that differentiation starts with large changes in H3K36me3 dynamics, and then relaxes towards the macrophage state.<\/p>\n

In summary, we applied scChIX-seq to two active histone modifications to find dynamic relationships between activation states. We then model these dynamics to infer chromatin velocity during macrophage differentiation.<\/p>\n

Discussion<\/h3>\n

Here, we demonstrate that scChIX-seq can deconvolve multiplexed histone modifications, expanding the number of histone marks that can be profiled in single cells. Using simulations, purified cell types and whole bone marrow, we demonstrate that scChIX-seq can accurately map several histone marks. To show how scChIX-seq can reveal unique biological insights in more challenging systems, we applied scChIX-seq to study H3K36me3 and H3K9me3 dynamics during mouse organogenesis to reveal the joint transcriptional and heterochromatin relationships in single cells. scChIX-seq can identify complex cell-type relationships between histone modifications, such as when several cell types can share a similar heterochromatin landscape. Finally, we applied scChIX-seq to two active marks during macrophage in vitro differentiation to quantify the relationship between two correlating marks. Importantly, scChIX-seq is flexible in which histone modifications can be used. The correlation structure between modifications is inferred from the model and therefore does not require a priori assumptions of specific features of the two modifications. Thus, scChIX-seq complements a recent method that focuses on differences in fragment lengths between Pol2 serine-5 phosphate and H3K27me3 to assign reads to their respective mark46<\/a><\/sup>.<\/p>\n

Recently, there have been new experimental innovations to CUT&TAG that modify the pA-Tn5 complex to map several histone modifications in single cells24<\/a>,47<\/a>,48<\/a>,49<\/a><\/sup>. One drawback of Tn5-based approaches (for example, CUT&TAG) compared with MNase-based (for example, sortChIC and CUT&RUN) used in this study is that Tn5 can have biases to open chromatin50<\/a><\/sup>. Current CUT&TAG methods suppress this bias by using more stringent washing conditions51<\/a><\/sup>, but exceedingly high salt conditions reduce the sensitivity and could wash away weakly bound factors such as transcription factors50<\/a>,51<\/a><\/sup>. On the flip side, MNase-based approaches involve more experimental effort than Tn5-based approaches, reducing the number of single cells that can be processed per round. Although we demonstrate our scChIX-seq method using an MNase-based approach (sortChIC), our computational and experimental framework can also be applied to Tn5-based strategies. Furthermore, our scChIX-seq method may have synergies with recent nanobody-based methods47<\/a>,48<\/a><\/sup>. For example, using two nanobodies, each specific to a different species of immunoglobulin G, one can profile four histone modifications by generating two sets of scChIX-seq simultaneously: two antibodies raised from one species and the other two antibodies raised from the second species.<\/p>\n

A limitation in scChIX-seq is that the maximum number of cuts at a specific base pair location is fundamentally limited by the copy number in that cell. Therefore, a nucleosome that has several modifications in their histone tails would still be cut only once. Currently, our binning strategy (5\u2009kilobase (kb), 50\u2009kb or gene bodies, depending on the biological question) and multinomial model assumes that there is no limit to the number of fragments that can be generated in one bin, which is an approximation that is valid when the bins are large and the number of cuts within the bins are small (for example, due to dropouts).<\/p>\n

We demonstrate that scChIX-seq can reveal biological insights by multimodal analysis that would otherwise be obscured by analyzing each modality separately. Overall, scChIX-seq unlocks multimodal analysis in antibody-based chromatin profiling and enables joint analysis of distinct histone modifications in single cells.<\/p>\n<\/div>\n<\/div>\n

\n

Methods<\/h2>\n
\n

Animal experiments<\/h3>\n

All mice used in this study were Cast-EiJ\/Bl6 mice and were bred and maintained in the Hubrecht Institute Animal Facility. All mouse experimentation was approved by the Animal Experimentation Committee (DEC) from the Koninklijke Nederlandse Akademie van Wetenschappen (KNAW) and complied with existing European Union legislation and local standards.<\/p>\n

Mouse bone marrow<\/h4>\n

Male 13-week-old C57BL\/6 mice were used to extract bone marrow cells. Femurs and tibia were extracted, the bone ends were cut away to access the bone marrow, which was flushed out using a 22G syringe with HBSS (\u2013\u2009calcium, \u2013\u2009magnesium, \u2013\u2009phenol red, Gibco, catalog no. 14175053) supplemented with Pen-Strep and 1% fetal calf serum. The bone marrow was dissociated and debris removed by passing through a 70\u2009\u03bcm cell strainer (Corning, catalog no. 431751). Cells were washed with 25\u2009ml supplemented HBSS before depleting the sample of unnucleated cells using IOTest 3 Lysing solution (Beckman Coulter) following the provider\u02bcs instructions. Cells were washed an additional two times with PBS before processing them by the sortChIC protocol for histone modifications. For whole bone marrow experiments (that is, not enriched for specific cell types), we processed cells using the sortChIC protocol for unfixed cells (without ethanol fixation). For the ground truth experiment with sorted cell types, we processed cells using the sortChIC protocol for ethanol-fixed cells. For ethanol fixation, cells were resuspended in 70% ethanol and fixed for 1\u2009h at \u201320\u2009\u00b0C. Afterwards cells were resuspended in Storage buffer (42.5\u2009ml H2<\/sub>O RNAse free, 1\u2009ml 1\u2009M HEPES pH\u20097.5 (Invitrogen), 1.5\u2009ml 5\u2009M NaCl, 3.6\u2009\u03bcl spermidine (Sigma Aldrich, catalog no. S2626-5G), protease inhibitor (Sigma Aldrich, catalog no. 5056489001), 200\u2009\u03bcl 0.5 M EDTA, 5\u2009\u03bcl dimethylsulfoxide) and frozen at \u201380\u2218<\/span><\/sup>C, before processing by the sortChIC protocol.<\/p>\n

Mouse organogenesis<\/h4>\n

No randomization or blinding was performed. Sex of embryos was not known at the time of collection. Four to five embryos were pooled for each reported timepoint (E9.5, E10.5, E11.5) before single-cell isolation. Pooled embryos were dissociated in TrypleE for 10\u2009min at room temperature. Undigested portions were physically removed and the remainder filtered through a 30\u2009\u03bcm filter before the single-cell suspension was split into three samples for each timepoint and each scChIX-seq experiment. Per timepoint, two single-cell samples were used each for a single antibody incubation (H3K36me3 or H3K9me3) and one sample for the double antibody incubation (H3K36me3 + H3K9me3). Antibody incubation was performed as described in the scChIX-seq protocol before single-cell capture using flow cytometry. A DNA library was prepared for each sample using the sortChIC protocol for unfixed cells.<\/p>\n

In vitro macrophage differentiation<\/h4>\n

For in vitro differentiation of bone marrow-derived macrophages, bone marrow was collected aseptically by flushing tibia and femurs from euthanized wild-type male C57BL\/6 mice with sterile RPMI and 10% FCS through a 70\u2009\u03bcm cell strainer (Corning). To enrich for stem and progenitor cells, lineage marker-positive (Lin+<\/sup>) cells were depleted by magnetic-activated cell sorting using a mouse Lineage Cell Depletion kit (Miltenyi Biotec). Lin\u2013<\/sup> cells were cultured on nontissue-culture-treated plates (Corning) for 7\u2009days in RPMI medium supplemented with 10% FCS, 100\u2009IU\u2009ml\u20131<\/sup> penicillin, 100\u2009mg\u2009ml\u20131<\/sup> streptomycin and 10\u2009ng\u2009ml\u20131<\/sup> recombinant murine MCSF (Peprotech). Medium was refreshed after 3\u2009days. Every 24\u2009h, suspension cells were collected and adherent cells were harvested by incubating 10\u2009min in 2\u2009mM EDTA\/0.5% BSA in PBS. Suspension and adherent cells were combined and stained with CellTrace fluorescent labels (Thermo Fisher), according to manufacturer\u2019s instructions. Briefly, cells were pelleted and resuspended in 37\u2009\u00b0C PBS containing fluorescent dyes (working concentrations CellTrace CSFE (CTC): 2.5\u2009\u03bcM; CellTrace Yellow (CTY): 2.5\u2009\u03bcM; CellTrace Far Red (CTFR): 0.5\u2009\u03bcM) at a concentration of 1,000,000 cells ml\u20131<\/sup>. Cells were incubated at 37\u2009\u00b0C protected from light for 20\u2009min. Staining reactions were stopped by adding two volumes of RPMI\/10% FCS and incubating for 5\u2009min at room temperature, protected from light, after which cells were washed twice in PBS. The following combinations of labels were used: unstained (day\u20090), CTC (day\u20091), CTY (day\u20092), CTFR (day\u20093), CTC\u2009+\u2009CTY (day\u20094), CTC\u2009+\u2009CTFR (day\u20095), CTY\u2009+\u2009CTFR (day\u20096) and CTC\u2009+\u2009CTY\u2009+\u2009CTFR (day\u20097). After harvesting and staining, cells were fixed in 70% ethanol for 1\u2009h and stored for later by the sortChIC protocol for fixed cells.<\/p>\n

Cell preparation without ethanol fixation for sortChIC experiments<\/h3>\n

Cells from whole bone marrow (H3K4me1+H3K27me3) and mouse embryos (H3K36me3+H3K9me3) were processed using the sortChIC method without ethanol fixation. Cells were processed in 0.5\u2009ml protein low-binding tubes. Following steps were performed on ice. Cells were resuspended in 500\u2009\u03bcl wash buffer (47.5\u2009ml H2<\/sub>O RNAse free, 1\u2009ml 1\u2009M HEPES pH\u20097.5 (Invitrogen), 1.5\u2009ml 5M NaCl, 3.6\u2009\u03bcl pure spermidine solution (Sigma Aldrich)). Cells were pelleted at 600g<\/i> for 3\u2009min and resuspended in 400\u2009\u03bcl wash buffer 1 (wash buffer with 0.05% saponin (Sigma Aldrich), protease inhibitor cocktail (Sigma Aldrich), 4\u2009\u03bcl 0.5\u2009M EDTA) containing the primary antibody (1:100 dilution for the antibody, saponin has to be prepared fresh every time as a 10% solution in PBS). Cells were incubated overnight at 4\u2009\u00b0C on a roller, before being washed once with 500\u2009\u03bcl wash buffer 2 (wash buffer with 0.05% saponin, protease inhibitor). Afterwards cells were resuspended in 500\u2009\u03bcl wash buffer 2 containing Protein A-Micrococcal Nuclease (pA-MNase) (3\u2009ng\u2009ml\u20131<\/sup>) and incubated for 1\u2009h at 4\u2009\u00b0C on a roller.<\/p>\n

Finally, cells were washed an additional two times with 500\u2009\u03bcl wash buffer 2 before passing it through a 70\u2009\u03bcm cell strainer (Corning, catalog no. 431751) and sorting G1 cells based on Hoechst staining on a BD Influx FACS machine into 384-well plates containing 50\u2009nl wash buffer 3 (wash buffer containing 0.05% saponin) and 5\u2009\u03bcl sterile filtered mineral oil (Sigma Aldrich) per well. Small volumes were distributed using a Nanodrop II system (Innovadyme).<\/p>\n

Cell preparation with ethanol fixation and surface antibody incubation for sortChIC experiments<\/h3>\n

Cells from sorted bone marrow (H3K27me3+H3K9me3) and macrophage in vitro differentiation (H3K4me1+H3K36me3) were processed using the ethanol fixation protocol. Sorted bone marrow cells were also incubated with surface antibody to enrich for known cell types. For the ethanol-fixed cells the above described sortChIC protocol was adapted. Wash buffers were used as described above, except that 0.05% saponin was exchanged for 0.05% Tween. Ethanol-fixed cells were thawed on ice. Cells were spun at 400g<\/i> for 5\u2009min and washed once with 400\u2009\u03bcl wash buffer 1. Cells were spun again at 400g<\/i> and resuspended in 400\u2009\u03bcl wash buffer 1. Cell suspension was split into three samples each having a volume of 400\u2009\u03bcl and incubated with one or two antibodies (1:100 dilution for the antibody) overnight on a roller at 4\u2009\u00b0C. The next day, cells were spun at 400g<\/i>, washed once with 400\u2009\u03bcl wash buffer 2 and resuspended in 500\u2009\u03bcl wash buffer 2 containing pA-MNase (3\u2009ng\u2009ml\u20131<\/sup>) and incubated for 1\u2009h on a rotator at 4\u2009\u00b0C. Next, cells were spun at 400g<\/i> and resuspended in 400\u2009\u03bcl wash buffer 2 (with addition of 5% blocking rat serum). To sort for defined cell types in the ground truth bone marrow experiment, surface antibodies were added according to these concentrations and were incubated for 30\u2009min on ice:<\/p>\n

\n

$$begin{array}{l}begin{array}{ll}{{mbox{antibody}}},&,{{mbox{info}}}\\ {{mbox{GR1}}},&,{{mbox{A647, anti-mouse Ly-6G\/Ly-6C (Gr-1) Antibody,}}}\\ & {mbox{clone: RB6-8C5}}\\ {{mbox{NK1}}},&,{{mbox{A488, anti-mouse NK-1.1 Antibody, clone: PK136}}}\\ {{mbox{CD19}}},&,{{mbox{BV421, anti-mouse CD19 Antibody, clone: 6D5}}}end{array}\\begin{array}{l}{{mbox{working concentration}}}\\1:8,000\\1:400\\1:200end{array}end{array}$$<\/span><\/p>\n<\/div>\n

BD FAC software v.1.2.0.142 was used to collect data from the FACS machine during cell sorting; see Supplemental Fig. 1<\/a> for the gating strategy.<\/p>\n

Finally, samples were washed once with 500\u2009\u03bcl wash buffer 2 before passing them through a 70\u2009\u03bcm cell strainer (Corning, catalog no. 431751) and sorting on a BD Influx FACS machine, with surface antibody specific gating, into 384-well plates containing 50\u2009nl wash buffer 3 (wash buffer containing 0.05% Tween) and 5\u2009\u03bcl sterile filtered mineral oil (Sigma Aldrich) per well. Small volumes were distributed using a Nanodrop II system (Innovadyme).<\/p>\n

MNase activation for sortChIC experiments<\/h3>\n

Targeted fragmentation was started by the addition of 5\u2009\u03bcl wash buffer 2 containing 4 mM CaCl2<\/sub>. For digestion, plates were incubated for 30\u2009min in a PCR machine set at 4\u2009\u00b0C. Afterwards the reaction was stopped by adding 100\u2009nl of a stop solution containing 40 mM EGTA, 1.5% NP40, and 10\u2009nl 2\u2009mg\u2009ml\u22121<\/sup> proteinase K. Plates were incubated in a PCR machine for further 20\u2009min at 4\u2009\u00b0C, before chromatin was released and pA-MNase permanently destroyed by proteinase K digestion at 65\u2009\u00b0C for 6\u2009h followed by 80\u2009\u00b0C for 20\u2009min to heat inactivate proteinase K. Afterwards plates were stored at \u201380\u2009\u00b0C until further processing.<\/p>\n

Library preparation for sortChIC experiments<\/h3>\n

DNA fragments were blunt ended by adding 150\u2009nl end repair mix per well and incubating for 30\u2009min at 37\u2009\u00b0C followed by 20\u2009min at 75\u2009\u00b0C for enzyme inactivation. End repair mix per well: Klenow large (NEB, catalog no. M0210L) 2.5\u2009nl, T4 PNK (NEB, catalog no. M0201L) 2.5\u2009nl, dNTPs 10\u2009mM 6\u2009nl, ATP 100\u2009mM 3.5\u2009nl, MgCl2<\/sub> 25\u2009mM 10\u2009nl, PEG8000 50% 7.5\u2009nl, PNK buffer 10\u00d7 (NEB, catalog no. B0201S) 35\u2009nl, BSA 20\u2009ng 1.8\u2009nl, nuclease-free water 81.3\u2009nl.<\/p>\n

Blunt fragments were subsequently A-tailed by adding 150\u2009nl per well of A-tailing mix and incubated for 15\u2009min at 72\u2009\u00b0C. Through the strong preference of AmpliTaq 360 to incorporate dATP as a single base overhang even in the presence of other nucleotides, a general dNTP removal was not necessary. A-tailing mix per well: AmpliTaq 360 (Thermo Fisher Scientific, catalog no. 4398828) 1\u2009nl, dATPs 100\u2009mM 1\u2009nl, KCl 1\u2009M 25\u2009nl, PEG8000 50% 7.5\u2009nl, BSA 20\u2009ng 0.8\u2009nl, nuclease-free water 114.8\u2009nl.<\/p>\n

Fragments were ligated to T-tail containing forked adapters containing a T7 polymerase binding site for in vitro transcription (IVT)-based amplification.<\/p>\n

Top strand: 5\u2032-GGTGATGCCGGTAATACGACTCACTATAGGGAGTTCTACAGTCCGACGATCNNNACACACTAT-3\u2032<\/p>\n

Bottom strand: 5\u2032-TAGTGTGTNNNGATCGTCGGACTGTAGAACTCCCTATAGTGAGTCGTATTACCGGCGAGCTT-3\u2032<\/p>\n

The three random nucleotides (NNN) were the unique molecular identifier used for read deduplication and the eight bases afterwards represent the cell barcodes, which were different for each of the 384 wells. For a full list of adapters and the cell barcodes for each well, see the excel sheet in Supplemental Table 1<\/a>. Cell barcodes for each 384-well plates are also found as a text file in the scChIX-seq Github repository: (https:\/\/github.com\/jakeyeung\/scChIX\/blob\/main\/inst\/extdata\/cellbarcodes_384_NLA_annotated.bc<\/a>).<\/p>\n

For ligation, 50\u2009nl of 5\u2009\u03bcM adapter in 50\u2009mM Tris pH\u20097 was added to each well with a Mosquito HTS (ttp labtech). After centrifugation, 150\u2009nl of ligation mix was added before incubating plates for 20\u2009min at 4\u2009\u00b0C, followed by 16\u2009h at 16\u2009\u00b0C for ligation and 10\u2009min at 65\u2009\u00b0C to inactivate ligase. Adapter ligation mix per well: T4 ligase (400,000\u2009U\u2009ml\u20131<\/sup>, NEB, catalog no. M0202L) 25\u2009nl, MgCl2<\/sub> 1\u2009M 3.5\u2009nl, Tris 1\u2009M pH\u20097.5 10.5\u2009nl, DTT 0.1\u2009M 52.5\u2009nl, ATP 100\u2009mM 3.5\u2009nl, PEG8000 50% 10\u2009nl, BSA 20\u2009ng 1\u2009nl, nuclease-free water 44\u2009nl.<\/p>\n

Before pooling, 1\u2009\u03bcl nuclease-free water was added to each well to minimize material loss. Ligation products were pooled by centrifugation into oil coated VBLOK200 Reservoir (ClickBio) at 500g<\/i> for 2\u2009min and the liquid face was transferred into 1.5\u2009ml Eppendorf tubes and then purified by centrifugation at 13,000g<\/i> for 1\u2009min and transferred into a fresh tube twice. DNA fragments were purified using Ampure XP beads (Beckman Coulter, prediluted one in eight in bead binding buffer: 1\u2009M NaCl, 20% PEG8000, 20\u2009mM Tris pH\u20098, 1\u2009mM EDTA) at a bead to sample ratio of 0.8. After 15\u2009min incubation at room temperature, beads were washed twice with 1\u2009ml 80% ethanol resuspending the beads during the first wash and resuspended in 20\u2009\u03bcl nuclease-free water. After 2\u2009min elution, the supernatant was transferred into a fresh 0.5\u2009ml tube. A second cleanup was performed adding 26\u2009\u03bcl undiluted Ampure XP beads and the beads were resuspended in 8\u2009\u03bcl nuclease-free water. The cleaned DNA was then linear amplified by IVT by adding 12\u2009\u03bcl of MEGAscript T7 Transcription Kit (Fisher Scientific, catalog no. AMB13345) for 12\u2009h at 37\u2009\u00b0C. Template DNA was removed by addition of 2\u2009\u03bcl\u20131<\/sup> TurboDNAse (IVT kit) and incubation for 15\u2009min at 37\u2009\u00b0C. The RNA produced was further purified using RNA Clean XP beads (Beckman Coulter) at a beads to sample ratio of 0.8 and samples were resuspended in 22\u2009\u03bcl of nuclease-free water. RNA was fragmented by mixing in 4.4\u2009\u03bcl fragmentation buffer (200 mM Tris-acetate pH\u20098.1, 500\u2009mM KOAc, 150\u2009mM MgOAc) and incubation for 2\u2009min at 94\u2009\u00b0C. Fragmentation was stopped by transferring samples to ice, adding 2.64\u2009\u03bcl 0.5\u2009M EDTA and another bead cleanup; samples were resuspended in 12\u2009\u03bcl nuclease-free water.<\/p>\n

RNA (5\u2009\u03bcl) was primed for reverse transcription by adding 0.5\u2009\u03bcl 10\u2009mM dNTPs and 1\u2009\u03bcl 20\u2009mM randomhexamerRT primer (5\u2032-GCCTTGGCACCCGAGAATTCCANNNNNN-3\u2032) and hybridizing it by incubation at 65\u2009\u00b0C for 5\u2009min followed by direct cool down on ice. Reverse transcription was performed by further addition of 2\u2009\u03bcl first strand buffer (part of Invitrogen kit, catalog no. 18064014), 1\u2009\u03bcl 0.1\u2009M DTT, 0.5\u2009\u03bcl RNAseOUT (Invitrogen, catalog no. LS10777019) and 0.5\u2009\u03bcl SuperscriptII (Invitrogen, catalog no. 18064014) and incubating the mixture at 25\u2009\u00b0C for 10\u2009min followed by 1\u2009h at 42\u2009\u00b0C. Single-stranded DNA was purified through incubation with 0.5\u2009\u03bcl RNAseA (Thermo Fisher, catalog no. EN0531) and incubation for 30\u2009min at 37\u2009\u00b0C.<\/p>\n

A final PCR amplification to add the Illumina small RNA barcodes and handles was performed by adding 25\u2009\u03bcl of NEBNext Ultra II Q5 Master Mix (NEB, catalog no. M0492L), 11\u2009\u03bcl nuclease-free water and 2\u2009\u03bcl of 10\u2009\u03bcM RP1 and RPIx primers.<\/p>\n

PCR protocol for sortChIC experiments<\/h3>\n

Activation for 30\u2009s at 98\u2009\u00b0C, 8\u201312 cycles (depending on starting material), 10\u2009s at 98\u2009\u00b0C, 30\u2009s at 60\u2009\u00b0C, 30\u2009s at 72\u2009\u00b0C, final amplification 10\u2009min at 72\u2009\u00b0C.<\/p>\n

PCR products were cleaned by two consecutive DNA bead clean-ups with a bead to sample ratio of 0.8. Final product was eluted in 7\u2009\u03bcl nuclease-free water. The abundance and quality of the final library were assessed by QUBIT and bioanalyzer.<\/p>\n

Data processing<\/h3>\n

All DNA libraries were sequenced on a Illumina NextSeq500 with 2\u2009\u00d7\u200975\u2009bp. We ran the raw fastq<\/span> files through the Single-Cell MultiOmics (SCMO) workflow (github.com\/BuysDB\/SingleCellMultiOmics<\/a>52<\/a><\/sup>). The workflow comprises of six steps.<\/p>\n

(1) Demultiplex raw fastq<\/span> files using demux.py<\/span> (SCMO). (2) Trim fastq<\/span> files by removing adapters using cutadapt<\/span> (v.3.5). (3) Map trimmed fastq<\/span> files using bwa<\/span> (v.0.7.17-r1188). (4) Tag bam files with cell barcode information, using bamtagmultiome.py<\/span> (SCMO). (5) Generate count tables using bamToCountTable.py<\/span> (SCMO). (6) Run dimensionality reduction of count matrices using run_LDA_model.R<\/span>. See an example of the pipeline in the scChIX-seq Github repository53<\/a><\/sup>.<\/p>\n

Unmixing scChIX-seq signal<\/h3>\n

Single-cell epigenomics techniques (for example, sortChIC, CUT&RUN and CUT&TAG) generate a vector of counts indicating the number of cut fragments that map in each genomic region for each cell. We model the vector of counts from a double-incubated cell (overrightarrow{y})<\/span> as a linear combination of two multinomial distributions: one coming from a cluster c<\/i> of histone modification 1, parameterized by ({overrightarrow{p}}_{c})<\/span>, the other from another cluster d<\/i> of histone modification 2 ({overrightarrow{q}}_{d})<\/span>. The log-likelihood for a linear combination of two multinomials is:<\/p>\n

\n

$${{{{rm{L}}}}}_{(c,d)}=log (Pleft(overrightarrow{y}| {overrightarrow{p}}_{c},{overrightarrow{q}}_{d},wright))propto mathop{sum }limits_{g=1}^{G}{y}_{g}log left(w{p}_{c,g}+left(1-wright){q}_{d,g}right).$$<\/span><\/p>\n

\n (1)\n <\/p>\n<\/div>\n

(overrightarrow{y})<\/span> is the number of cuts across the genome for a double-incubated cell. p<\/i>c<\/i>,g<\/i><\/sub> and q<\/i>d<\/i>,g<\/i><\/sub> are cluster-specific probabilities indicating the likelihood that a cut fragment maps to region g<\/i> in histone modifications 1 and 2, respectively. w<\/i> is the mixing fraction of histone modification 1 in the double-incubated cell, which we estimate by maximizing the log-likelihood given (overrightarrow{y})<\/span>, ({overrightarrow{p}}_{c})<\/span> and ({overrightarrow{q}}_{d})<\/span>.<\/p>\n

Applying single-cell techniques to complex tissues generates data with many clusters. Therefore, given a double-incubated cell, we do not know which pair of clusters (c<\/i>,d<\/i>) were combined to generate the observed counts. We therefore calculate the log-likelihood for all possible pairs of clusters learned from the training data and then select the cluster pair with the highest probability for each cell.<\/p>\n

Cluster-specific probabilities ({overrightarrow{p}}_{c})<\/span> and ({overrightarrow{q}}_{d})<\/span> are learned by applying LDA (with k\u2009=\u200930 topics) using the topicmodels<\/span> R package54<\/a><\/sup> to the training data (that is, single-incubated cells), which are count matrices.<\/p>\n

After assigning each cell to the most probable cluster pair ((hat{c},hat{d}))<\/span>, we assign y<\/i>i,j<\/sub>, the jth read mapped to region g<\/i> in cell i<\/i>, to histone mark 1 with probability P<\/i>i,j<\/sub>:<\/p>\n

\n

$${P}_{mathrm{i,j}}=frac{w{p}_{hat{c},g}}{w{p}_{hat{c},g}+left(1-wright){q}_{hat{d},g}}.$$<\/span><\/p>\n

\n (2)\n <\/p>\n<\/div>\n

This assignment generates a pair of vectors ({overrightarrow{y}}_{1,i})<\/span> and ({overrightarrow{y}}_{2,i})<\/span> that are linked because they both come from cell i<\/i>. Unmixed counts ({overrightarrow{y}}_{1,i})<\/span> and ({overrightarrow{y}}_{2,i})<\/span> are then projected back onto the space inferred from training data of histone modification 1 and 2, respectively. The links between histone modification 1 and 2 are used to transfer labels and create linked UMAPs between the two histone modifications.<\/p>\n

Latent Dirichlet allocation<\/h3>\n

LDA is a probabilistic matrix decomposition model that is useful when the input data is a matrix of counts. LDA uses hierarchical multinomial models to estimate the relative frequencies of cuts in each genomic region in single cells.<\/p>\n

To generate the genomic location of the jth read for cell i:<\/p>\n

Choose a topic z<\/i>i,j<\/sub> by sampling from the cell-specific distribution of topics:<\/p>\n

\n

$$begin{array}{r}{overrightarrow{U}}_{mathrm{i}} sim ,{{{rm{Dirichlet}}}},(alpha )\\ {z}_{mathrm{i,j}} sim ,{{{rm{Multinomial}}}},({overrightarrow{U}}_{i},1)end{array}$$<\/span><\/p>\n<\/div>\n

Choose genomic region w<\/i>i,j<\/sub> by sampling from the topic-specific distribution of genomic regions:<\/p>\n

\n

$$begin{array}{r}{overrightarrow{V}}_{mathrm{k}} sim ,{{{rm{Dirichlet}}}},(delta )\\ {w}_{mathrm{i,j}} sim ,{{{rm{Multinomial}}}},({overrightarrow{V}}_{{z}_{mathrm{i,j}}},1)end{array}$$<\/span><\/p>\n<\/div>\n

The Dirichlet distributions are priors to prevent overfitting when there are few cuts in the region. We used the LDA model implemented by the topicmodels R package, using the Gibbs sampling implementation with hyperparameters \u03b1<\/i>\u2009=\u20091.67, \u03b4<\/i>\u2009=\u20090.1, where K is the number of topics23<\/a><\/sup>.<\/p>\n

We estimate ({overrightarrow{p}}_{c})<\/span> and ({overrightarrow{q}}_{d})<\/span> for each cluster in histone modification 1 ({{overrightarrow{p}}_{1},{overrightarrow{p}}_{2},…,{overrightarrow{p}}_{C}})<\/span> and modification 2 ({{overrightarrow{q}}_{1},{overrightarrow{q}}_{2},…,{overrightarrow{q}}_{D}})<\/span> by averaging the estimated probabilities across cells assigned to each cluster for each gene g<\/i>:<\/p>\n

\n

$$p_{g,c}=frac{1}{vert C vert}mathop{sum }limits_{mathrm{i}in C}mathop{sum }limits_{mathrm{k}=1}^{K}{V}_{mathrm{g,k}}{U}_{mathrm{k,i}}$$<\/span><\/p>\n<\/div>\n

where C<\/i> is the set of cells that belong to cluster c<\/i>.<\/p>\n

Simulation of single- and double-incubated histone modification data<\/h3>\n

To simulate multimodal single-cell histone modification data with varying degrees of overlap, we extended simATAC55<\/a><\/sup> to allow generating cell-type profiles from histone modifications of varying mutually exclusive relationships.<\/p>\n

For each cell type, we first run simATAC to generate sparse count data of 10,000 loci across 750 cells partitioned into three technical replicates of 250 cells each. The high-dimensional count data are sparse. Counts from each locus are generated according to a Poisson likelihood with locus-specific means (\u03bb<\/i>) matching real single-cell ATAC-seq from K562 cells (GSE99172).<\/p>\n

In our 750 cells, cells 1\u2013250 represent single-incubated cells from mark 1; cells 251\u2013500 from mark 2; cells 501\u2013750 from double-incubated cells. Cells from mark 1 have counts generated from locus-specific means \u03bb<\/i>. Cells from mark 2 also have counts generated from \u03bb<\/i>, but we swap the top x<\/i>% of bins with highest \u03bb<\/i> with bins with lowest \u03bb<\/i>, allowing precisely defined sets of mutually exclusive and overlapping bins. We use x<\/i>\u2009=\u20091%, 50% and 99% to benchmark our method from mostly overlapping (that is x<\/i>\u2009=\u20091%) to mostly mutually exclusive (that is x<\/i>\u2009=\u200999%) Cells from mark 3 are generated by adding counts generated from mark 1 and mark 2 to simulate double-incubated cells.<\/p>\n

To generate cell-type-specific profiles, we repeat the above with a cell-type-specific random seed and shuffle the order of the bins. This generates count data where \u03bb<\/i> is cell-type specific, but the distribution of \u03bb<\/i> are preserved genome-wide.<\/p>\n

Estimating the top cluster-specific bins<\/h3>\n

We use the LDA matrix factorization to identify the top cluster-specific bins in the data. We rank the bin loadings for each cell type and take the top 150 (whole bone marrow) or 250 (mouse organogenesis) bins with the largest loadings.<\/p>\n

Inferring pseudotime in differentiation data<\/h3>\n

To analyze the macrophage differentiation data, we first removed erythroblasts, plasmacytoid dendritic cells, and innate lymphocyte cells from the data, which were concentrated at day\u20090 and not considered to be part of the macrophage differentiation trajectory. We then ran LDA (k\u2009=\u200930 topics) and performed principal component analysis (PCA) on the LDA outputs, which retrieves the principal components that explain the largest amount of variance after denoising the data. We used the first principal component for H3K4me1 and H3K36me3 to define pseudotime, which we found correlates with the day along the timecourse.<\/p>\n

Unmixing scChIX-seq signal from continuous pseudotime<\/h3>\n

To apply scChIX-seq on continuous pseudotime, we modify the log-likelihood (equation (1<\/a>)) to account for a continuous variable:<\/p>\n

\n

$${{{rm{L}}}}left({t}_{1},{t}_{2}right)=log left(Pleft(overrightarrow{y}| overrightarrow{p}left({t}_{1}right),overrightarrow{q}left({t}_{2}right),wright)right)propto mathop{sum }limits_{g=1}^{G}{y}_{g}log left(w{p}_{g}left({t}_{1}right)+left(1-wright){q}_{g}left({t}_{2}right)right)$$<\/span><\/p>\n

\n (3)\n <\/p>\n<\/div>\n

where t<\/i>1<\/sub>\u2009\u2208<\/span>\u2009[0,\u20091] is pseudotime from histone modification 1 and t<\/i>2<\/sub>\u2009\u2208<\/span>\u2009[0,\u20091] is pseudotime from modification 2.<\/p>\n

To estimate pseudotime, we ran LDA to denoise the count matrix, and then ran PCA to estimate largest principal components explaining the variance in the data. We took the first principal component as our pseudotime estimate for both marks, which captured the epigenomic changes over the 7-day timecourse.<\/p>\n

({p}_{g}left(tright))<\/span> is estimated by fitting the signal from histone modification 1 at genomic region g<\/i> with a lowess curve along pseudotime. We estimate q<\/i>g<\/i><\/sub> analogously but using signal from histone modification 2.<\/p>\n

To infer the pseudotime of histone modifications 1 and 2 simultaneously given a vector of counts from a double-incubated cell, we estimate t<\/i>1<\/sub> and t<\/i>2<\/sub> that minimizes the log-likelihood L<\/i> from equation (3<\/a>). We estimate the variance-covariance matrix of t<\/i>1<\/sub> and t<\/i>2<\/sub> by the square root of the inverse of the Hessian matrix, which we use to calculate the standard errors.<\/p>\n

Since the t<\/i>1<\/sub> and t<\/i>2<\/sub> are constrained between 0 and 1, we use the L-BFGS-B optimization algorithm implemented in R. Since estimates from a single cell can sometimes be noisy due to low counts, we sum the counts across the 25-nearest neighbors (estimated from the latent space inferred by LDA) for each double-incubated cell.<\/p>\n

Chromatin velocity during macrophage differentiation<\/h3>\n

We assume that dynamic genomic regions in H3K36me3 can be modeled using a first-order differential equation<\/p>\n

\n

$$frac{d{K}_{36}left(tright)}{dt}={K}_{4}left(tright)-gamma {K}_{36}left(tright).$$<\/span><\/p>\n

\n (4)\n <\/p>\n<\/div>\n

We estimate the time constant \u03b3<\/i> for each genomic region by fitting an exponential relaxation function across pseudotime<\/p>\n

\n

$${K}_{36}left(tright)={y}_{0}+Aleft(1-{e}^{-gamma t}right),$$<\/span><\/p>\n

\n (5)\n <\/p>\n<\/div>\n

where y<\/i>0<\/sub> is the signal at t<\/i>\u2009=\u20090 and A<\/i> is the predicted H3K36me3 levels at steady state. Fitting the \u03b3<\/i> directly from the pseudotime allows us to leverage signal from both single- and deconvolved cells.<\/p>\n

To predict future values of H3K36me3 levels for each cell at each genomic region, we use the Euler method and plug in the estimated \u03b3<\/i>, H3K4me1 levels at time t<\/i> and time step h<\/i> of 0.02 pseudotime units:<\/p>\n

\n

$${K}_{36}left(t+1right)={K}_{36}left(tright)+hleft({K}_{4}left(tright)-gamma {K}_{36}left(tright)right).$$<\/span><\/p>\n

\n (6)\n <\/p>\n<\/div>\n

Finally, we project the single- and double-incubated H3K36me3 signal onto the first two principal components and project the predicted future values onto the PCA. We use the velocity grid flow visualization as implemented in velocyto56<\/a><\/sup> to visualize the velocity vectors on the PCA space.<\/p>\n

Comparison with multi-CUT&TAG<\/h3>\n

Raw fastq files (R1, R2 and R3) from the single-cell experiments were downloaded from Gene Expression Omnibus accession number GSE171554. The first 42 bases of the reads in R1 and R2 were trimmed to remove the barcodes and the bases common to all Tn5 adapter sequences. The 16-base cell barcodes in R3 were added to the fastq headers of R1 and R2. The trimmed and cell-barcoded R1 and R2 reads were then aligned to the mm10 mouse genome using Burrows-Wheeler aligner (bwa v.0.7.17-r1188). Fragments that start at same location and have the same cell barcode were considered duplicates and discarded. Cells with more than 100 fragments with MAPQ scores in R1 greater than or equal to 40 were kept for comparison with scChIX-seq.<\/p>\n

Reporting summary<\/h3>\n

Further information on research design is available in the Nature Portfolio Reporting Summary<\/a> linked to this article.<\/p>\n<\/div>\n<\/div><\/div>\n

\n
\n

Data availability<\/h2>\n

The data discussed in this publication have been deposited in NCBI\u2019s Gene Expression Omnibus and are accessible through Gene Expression Omnibus Series accession number GSE155280<\/a> (ref. 57<\/a><\/sup>).<\/p>\n<\/div>\n

\n

Code availability<\/h2>\n
\n

We developed the SingleCellMultiOmics package, in which there are modules used for processing sortChIC data (https:\/\/github.com\/BuysDB\/SingleCellMultiOmics)<\/a>52<\/a><\/sup>, and an R package that implements scChIX-seq and contains snakemake workflows for processing data and example notebooks for downstream analyses (https:\/\/github.com\/jakeyeung\/scChIX<\/a>)53<\/a><\/sup>.<\/p>\n<\/p><\/div>\n<\/div>\n

\n

References<\/h2>\n
\n
    \n
  1. \n

    Rothbart, S. B. & Strahl, B. D. Interpreting the language of histone and DNA modifications. Biochim. Biophys. Acta<\/i> 1839<\/b>, 627\u2013643 (2014).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  2. \n

    Kundaje, A. et al. Integrative analysis of 111 reference human epigenomes. Nature<\/i> 518<\/b>, 317\u2013330 (2015).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  3. \n

    Landt, S. G. et al. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res.<\/i> 22<\/b>, 1813\u20131831 (2012).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  4. \n

    Yue, F. et al. A comparative encyclopedia of DNA elements in the mouse genome. Nature<\/i> 515<\/b>, 355\u2013364 (2014).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  5. \n

    Lara-Astiaso, D. et al. Chromatin state dynamics during blood formation. Science<\/i> 345<\/b>, 943\u2013949 (2014).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  6. \n

    Rotem, A. et al. Single-cell ChIP\u2013seq reveals cell subpopulations defined by chromatin state. Nat. Biotechnol.<\/i> 33<\/b>, 1165\u20131172 (2015).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  7. \n

    Grosselin, K. et al. High-throughput single-cell ChIP\u2013seq identifies heterogeneity of chromatin states in breast cancer. Nat. Genet.<\/i> 51<\/b>, 1060\u20131066 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  8. \n

    Ai, S. et al. Profiling chromatin states using single-cell itChIP\u2013seq. Nat. Cell Biol.<\/i> 21<\/b>, 1164\u20131172 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  9. \n

    Schmid, M., Durussel, T. & Laemmli, U. K. ChIC and ChEC: genomic mapping of chromatin proteins. Mol. Cell<\/i> 16<\/b>, 147\u2013157 (2004).<\/p>\n

    CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  10. \n

    Skene, P. J., Henikoff, J. G. & Henikoff, S. Targeted in situ genome-wide profiling with high efficiency for low cell numbers. Nat. Protoc.<\/i> 13<\/b>, 1006\u20131019 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  11. \n

    Kaya-Okur, H. S. et al. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nat. Commun.<\/i> 10<\/b>, 1\u201310 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  12. \n

    Ku, W. L. et al. Single-cell chromatin immunocleavage sequencing (scChIC-seq) to profile histone modification. Nat. Methods<\/i> 16<\/b>, 323\u2013325 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  13. \n

    Wang, Q. et al. CoBATCH for high-throughput single-cell epigenomic profiling. Mol. Cell<\/i> 76<\/b>, 206\u2013216.e7 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  14. \n

    Harada, A. et al. A chromatin integration labelling method enables epigenomic profiling with lower input. Nat. Cell Biol.<\/i> 21<\/b>, 287\u2013296 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  15. \n

    Wu, S. J. et al. Single-cell CUT&Tag analysis of chromatin modifications in differentiation and tumor progression. Nat. Biotechnol.<\/i> 39<\/b>, 819\u2013824 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  16. \n

    Bartosovic, M., Kabbe, M. & Castelo-Branco, G. Single-cell CUT&Tag profiles histone modifications and transcription factors in complex tissues. Nat. Biotechnol.<\/i> 39<\/b>, 825\u2013835 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  17. \n

    Janssens, D. H. et al. Automated CUT&Tag profiling of chromatin heterogeneity in mixed-lineage leukemia. Nat. Genet.<\/i> 53<\/b>, 1586\u20131596 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  18. \n

    Zeller, P. et al. Hierarchical chromatin regulation during blood formation uncovered by single-cell sortChIC. Preprint at bioRxiv<\/i> https:\/\/doi.org\/10.1101\/2021.04.26.440606<\/a> (2021).<\/p>\n<\/li>\n

  19. \n

    Ku, W. L., Pan, L., Cao, Y., Gao, W. & Zhao, K. Profiling single-cell histone modifications using indexing chromatin immunocleavage sequencing. Genome Res.<\/i> 31<\/b>, 1831\u20131842 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  20. \n

    Navidi, Z., Zhang, L. & Wang, B. simATAC: a single-cell ATAC-seq simulation framework. Genome Biol.<\/i> 22<\/b>, 1\u201316 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  21. \n

    Pauler, F. M. et al. H3K27me3 forms BLOCs over silent genes and intergenic regions and specifies a histone banding pattern on a mouse autosomal chromosome. Genome Res.<\/i> 19<\/b>, 221\u2013233 (2009).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  22. \n

    Blei, D. M., Ng, A. Y. & Jordan, M. I. Latent Dirichlet allocation. J. Mach. Learn. Res.<\/i> 3<\/b>, 993\u20131022 (2003).<\/p>\n


    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  23. \n

    Gr\u00fcn, B. & Hornik, K. topicmodels: an R package for fitting topic models. J. Stat. Softw.<\/i> 40<\/b>, 1\u201330 (2011).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  24. \n

    Gopalan, S., Wang, Y., Harper, N. W., Garber, M. & Fazzio, T. G. Simultaneous profiling of multiple chromatin proteins in the same cells. Mol. Cell<\/i> 81<\/b>, 4736\u20134746.e5 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  25. \n

    Giladi, A. et al. Single-cell characterization of haematopoietic progenitors and their trajectories in homeostasis and perturbed haematopoiesis. Nat. Cell Biol.<\/i> 20<\/b>, 836\u2013846 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  26. \n

    Mendenhall, E. M. et al. GC-rich sequence elements recruit PRC2 in mammalian ES cells. PLoS Genet.<\/i> 6<\/b>, e1001244 (2010).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  27. \n

    Zou, F. et al. Expression and function of tetraspanins and their interacting partners in B cells. Front. Immunol.<\/i> 9<\/b>, 1606 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  28. \n

    Benhamou, D. et al. The c-Myc\/miR17-92\/PTEN axis tunes PI3K activity to control expression of recombination activating genes in early B cell development. Front. Immunol.<\/i> 9<\/b>, 2715 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  29. \n

    Goldmit, M. et al. Epigenetic ontogeny of the Igk locus during B cell development. Nat. Immunol.<\/i> 6<\/b>, 198\u2013203 (2005).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  30. \n

    Pan, C., Baumgarth, N. & Parnes, J. R. CD72-deficient mice reveal nonredundant roles of CD72 in B cell development and activation. Immunity<\/i> 11<\/b>, 495\u2013506 (1999).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  31. \n

    Gr\u00fcn, D. et al. De novo prediction of stem cell identity using single-cell transcriptome data. Cell Stem Cell<\/i> 19<\/b>, 266\u2013277 (2016).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  32. \n

    Pishesha, N. et al. Transcriptional divergence and conservation of human and mouse erythropoiesis. Proc. Natl Acad. Sci. USA<\/i> 111<\/b>, 4103\u20134108 (2014).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  33. \n

    Koretzky, G. A., Abtahian, F. & Silverman, M. A. SLP76 and SLP65: complex regulation of signalling in lymphocytes and beyond. Nat. Rev. Immunol.<\/i> 6<\/b>, 67\u201378 (2006).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  34. \n

    Brachtendorf, G. et al. Early expression of endomucin on endothelium of the mouse embryo and on putative hematopoietic clusters in the dorsal aorta. Dev. Dyn.<\/i> 222<\/b>, 410\u2013419 (2001).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  35. \n

    Sedykh, I. et al. Zebrafish Rfx4 controls dorsal and ventral midline formation in the neural tube. Dev. Dyn.<\/i> 247<\/b>, 650\u2013659 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  36. \n

    DeBoer, E. M. et al. Prenatal deletion of the RNA-binding protein HuD disrupts postnatal cortical circuit maturation and behavior. J. Neurosci.<\/i> 34<\/b>, 3674\u20133686 (2014).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  37. \n

    Inoue, T. et al. Analysis of mouse Cdh6 gene regulation by transgenesis of modified bacterial artificial chromosomes. Dev. Biol.<\/i> 315<\/b>, 506\u2013520 (2008).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  38. \n

    Chen, A. F. et al. GRHL2-dependent enhancer switching maintains a pluripotent stem cell transcriptional subnetwork after exit from naive pluripotency. Cell Stem Cell<\/i> 23<\/b>, 226\u2013238.e4 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  39. \n

    Logan, M. et al. Expression of Cre recombinase in the developing mouse limb bud driven by aPrxl enhancer. Genesis<\/i> 33<\/b>, 77\u201380 (2002).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  40. \n

    Takeuchi, J. K. & Bruneau, B. G. Directed transdifferentiation of mouse mesoderm to heart tissue by defined factors. Nature<\/i> 459<\/b>, 708\u2013711 (2009).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  41. \n

    Zhao, R. et al. Loss of both GATA4 and GATA6 blocks cardiac myocyte differentiation and results in acardia in mice. Dev. Biol.<\/i> 317<\/b>, 614\u2013619 (2008).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  42. \n

    Cao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature<\/i> 566<\/b>, 496\u2013502 (2019).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  43. \n

    Murphy, Z. C. et al. Regulation of RNA polymerase II activity is essential for terminal erythroid maturation. Blood<\/i> 138<\/b>, 1740\u20131756 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  44. \n

    Gautier, E. L. et al. Gene-expression profiles and transcriptional regulatory pathways that underlie the identity and diversity of mouse tissue macrophages. Nat. Immunol.<\/i> 13<\/b>, 1118\u20131128 (2012).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  45. \n

    La Manno, G. et al. RNA velocity of single cells. Nature<\/i> 560<\/b>, 494\u2013498 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  46. \n

    Janssens, D. H. et al. CUT&Tag2for1: a modified method for simultaneous profiling of the accessible and silenced regulome in single cells. Genome Biol.<\/i> 23<\/b>, 81 (2022).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  47. \n

    Stuart, T. et al. Nanobody-tethered transposition allows for multifactorial chromatin profiling at single-cell resolution. Preprint at bioRxiv<\/i> https:\/\/doi.org\/10.1101\/2022.03.08.483436<\/a> (2022).<\/p>\n<\/li>\n

  48. \n

    Bartosovic, M. & Castelo-Branco, G. Multimodal chromatin profiling using nanobody-based single-cell CUT&Tag. Preprint at bioRxiv<\/i> https:\/\/doi.org\/10.1101\/2022.03.08.483459<\/a> (2022).<\/p>\n<\/li>\n

  49. \n

    Meers, M. P., Llagas, G., Janssens, D. H., Codomo, C. A. & Henikoff, S. Multifactorial profiling of epigenetic landscapes at single-cell resolution using MulTI-Tag. Nat. Biotechnol.<\/i> https:\/\/doi.org\/10.1038\/s41587-022-01522-9<\/a> (2022).<\/p>\n<\/li>\n

  50. \n

    Wang, M. & Zhang, Y. Tn5 transposase-based epigenomic profiling methods are prone to open chromatin bias. Preprint at bioRxiv<\/i> https:\/\/doi.org\/10.1101\/2021.07.09.451758<\/a> (2021).<\/p>\n<\/li>\n

  51. \n

    Kaya-Okur, H. S., Janssens, D. H., Henikoff, J. G., Ahmad, K. & Henikoff, S. Efficient low-cost chromatin profiling with CUT&Tag. Nat. Protoc.<\/i> 15<\/b>, 3264\u20133283 (2020).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  52. \n

    de Barbanson, B. A. et al. BuysDB\/SingleCellMultiOmics: 0.1.30 (v.0.1.30). Zenodo. https:\/\/doi.org\/10.5281\/zenodo.7074511<\/a> (2022).<\/p>\n<\/li>\n

  53. \n

    Yeung, J. jakeyeung\/scChIX: v.1.0.1 (v.1.0.1). Zenodo. https:\/\/doi.org\/10.5281\/zenodo.7152037<\/a> (2022).<\/p>\n<\/li>\n

  54. \n

    Gr\u00fcn, B. & Hornik, K. topicmodels: an R package for fitting topic models. J. Stat. Softw.<\/i> 40<\/b>, 1\u201330 (2011).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  55. \n

    Navidi, Z., Zhang, L. & Wang, B. simATAC: a single-cell ATAC-seq simulation framework. Genome Biol.<\/i> 22<\/b>, 74 (2021).<\/p>\n

    Article<\/a>\u00a0
    \n     CAS<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  56. \n

    La Manno, G. et al. RNA velocity of single cells. Nature<\/i> 560<\/b>, 494\u2013498 (2018).<\/p>\n

    Article<\/a>\u00a0
    \n
    \n Google Scholar<\/a>\u00a0\n <\/p>\n<\/li>\n

  57. \n

    Yeung, J., Florescu, M., Zeller, P, de Barbanson, B. A., Wellenstein, M. D. & van Oudenaarden, A. scChIX-seq infers relationships between histone modifications in single cells. Datasets. Gene Expression Omnibus. https:\/\/www.ncbi.nlm.nih.gov\/geo\/query\/acc.cgi?acc=GSE155280<\/a> (2022).<\/p>\n<\/li>\n<\/ol>\n

    Download references<\/a><\/p>\n<\/div>\n<\/div>\n

    \n

    Acknowledgements<\/h2>\n

    We thank M. van Loenhout for experimental advice on purifying cell types from the bone marrow, R. van der Linden for expertise with FACS and M. Blotenburg for help with cell typing the mouse organogenesis dataset. We thank M. Saraswat and O. Stegle for discussions on multinomial distributions. This work was supported by a European Research Council Advanced grant (ERC-AdG 742225-IntScOmics); Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) TOP grant (NWO CW 714.016.001) and NWO grant (OCENW.GROOT.2019.017); the Swiss National Science Foundation Early Postdoc Mobility (P2ELP3-184488 to P.Z. and P2BSP3-174991 to J.Y.); Marie Sklodowska-Curie Actions Postdoc (798573 to P.Z.) and the Human Frontier for Science Program Long-Term Fellowships (LT000209-2018-L to P.Z. and LT000097-2019-L to J.Y.). This work is part of the Oncode Institute which is financed partly by the Dutch Cancer Society.<\/p>\n<\/div>\n

    \n

    Author information<\/h2>\n
    \n

    Author notes<\/span><\/p>\n

      \n
    1. \n

      These authors contributed equally: Jake Yeung, Maria Florescu, Peter Zeller.<\/p>\n<\/li>\n<\/ol>\n

      Authors and Affiliations<\/h3>\n
        \n
      1. \n

        Oncode Institute, Hubrecht Institute-KNAW (Royal Netherlands Academy of Arts and Sciences) and University Medical Center Utrecht, Utrecht, the Netherlands<\/p>\n

        Jake Yeung,\u00a0Maria Florescu,\u00a0Peter Zeller,\u00a0Buys Anton de Barbanson,\u00a0Max D. Wellenstein\u00a0&\u00a0Alexander van Oudenaarden<\/p>\n<\/li>\n

      2. \n

        Institute of Science and Technology Austria (ISTA), Klosterneuburg, Austria<\/p>\n

        Jake Yeung<\/p>\n<\/li>\n<\/ol>\n

        Contributions<\/h3>\n

        J.Y., M.F., B.A.d.B. and A.v.O. conceived the project. M.F. developed double-incubation techniques and performed mouse bone marrow and organogenesis experiments with help from P.Z. P.Z. developed single-incubation techniques. P.Z. and M.D.W. designed and performed macrophage in vitro differentiation experiments. J.Y., M.F. and A.v.O. analyzed the data. J.Y. developed and applied statistical methods with help from M.F. and B.A.d.B. B.A.d.B. wrote the sortChIC preprocessing pipeline, with help from M.F. and J.Y. J.Y., M.F. and A.v.O. wrote the manuscript, with input from P.Z., M.D.W. and B.A.d.B.<\/p>\n

        Corresponding authors<\/h3>\n

        Correspondence to
        \n Jake Yeung<\/a> or Alexander van Oudenaarden<\/a>.<\/p>\n<\/div>\n<\/div>\n

        \n

        Ethics declarations<\/h2>\n
        \n

        Competing interests<\/h3>\n

        The authors declare no competing interests.<\/p>\n<\/p><\/div>\n<\/div>\n

        \n

        Additional information<\/h2>\n

        Publisher\u2019s note<\/b> Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.<\/p>\n<\/div>\n

        \n

        Extended data<\/h2>\n
        \n
        \n

        Extended Data Fig. 1 Benchmarking scChIX-seq across a range of overlapping patterns.<\/a><\/h3>\n

        Left column: simulation results in a mutually exclusive scenario (that is 1% of loci are overlapping). Middle column: results for an intermediate amount of overlap (that is 50% of loci are overlapping). Right column: results for highly correlated scenario (that is 99% of loci are overlapping). (a)<\/b> Distribution of unique fragment cuts per cell in simulation. (b)<\/b> Sparsity of the input matrix. Note that in the mutually exclusive scenario, the double-incubated marks is less sparse than single-incubated marks because loci with zero reads in one mark often have non-zero reads in another mark. (c)<\/b> Distribution of the degree of overlap (defined as fraction of double-incubated signal belonging to mark1: (p=frac{{S}_{1}}{{S}_{1}+{S}_{2}})<\/span>) for each locus genome-wide. (d)<\/b> Estimated degree of overlap from scChIX-seq. (e)<\/b> UMAP representation of the three cell types underlying simulation. UMAPs from the two marks are linked by double-incubated cells that are deconvolved by scChIX-seq. (f)<\/b> Empirical 95% confidence interval across the range of (hat{p}=frac{{hat{S}}_{1}}{{hat{S}}_{1}+{hat{S}}_{2}})<\/span> (from 0 to 1). Range obtained by aggregating results from the three overlapping patterns. n=101 simulation datapoints spread evenly between 0 and 1 inclusive. Error bars are empirial 95% confidence intervals, centers are the mean.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 2 scChIX-seq accurately deconvolves double-incubated signal into their respective histone modifications.<\/a><\/h3>\n

        (a)<\/b> Histogram of unique fragment cuts per cell. (b)<\/b> Histogram of fraction of unique fragments starting with a “TA” motif. (c, d)<\/b> UMAP of latent Dirichlet allocation (LDA) embedding using k=30 topics for H3K27me3 (c) and H3K9me3 (d). (e, f)<\/b> UMAP representation of H3K27me3 (left) and H3K9me3 (right) data colored by unmixed or single-incubated cells (e) or ground truth cell type labels defined by FACS (f). (g, h)<\/b> Genome-wide Pearson correlation between deconvolved H3K27me3 (g) and H3K9me3 (h) signal versus ground truth sortChIC purified by FACS. Shared genomic regions were calculated by using 1 kb bins across the genome. (i)<\/b> Comparison of fragments per cell obtained from Multi-CUT&TAG versus scChIX-seq. Multi-CUT&TAG data came from a mixture of embryonic and trophoblast stem cells in vitro<\/i>, while scChIX-seq came from sorted bone marrow cells in vivo<\/i>. n=1806 cells for Multi-CUT&TAG, n=290 for scChIX-seq.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 3 Coverage tracks of deconvolved cells and genome statistics.<\/a><\/h3>\n

        (a)<\/b> Coverage tracks for B cells visualizing the H3K27me3+H3K9me3, deconvolved H3K27me3 or H3K9me3, and ground truth H3K27me3 or H3K9me3 histone modification levels for three different genomic regions. Double-incubated signal in grey, H3K27me3 single, and unmixed signal in orange, and H3K9me3 single and unmixed signal in blue. Under each coverage track are cut fragments of single cells. Each row of the single cells track are cuts from an individual cell. Shown are a subset of cells, which were chosen for their high number of cuts in the region. Rows are in decreasing order of total number of cuts. (b)<\/b> H3K27me3 coverage tracks showing the region around Pax5<\/i> for the ground truth H3K27me3 pseudobulk signal from single-incubated cells and for the deconvolved H3K27me3 pseudobulk signal from double-incubated cells for three cell types: B cells (grey), granulocytes (green), and NK cells (blue). (c)<\/b> H3K9me3 (top) and H3K27me3 (bottom) coverage tracks showing the region around Auts2<\/i> for ground truth (single-incubated) and for the unmixed (unmixed) for B cells (grey), granulocytes (green) and NK cells (blue), respectively. d<\/b> Distribution of assignment probability estimates in the genome for the three cell types. Vertical dotted lines represent cutoffs to define H3K9me3-specific (that is p<\/i>\u2009<\u20090.5) or H3K27me3-specific regions (that is p<\/i>\u22650.5). e<\/b> Boxplot distributions of GC content in H3K27me3-marked and H3K9me3-marked regions. f<\/b> Boxplot distributions of distance to TSS in the two classes of regions. Distances are measured from the center of the 50 kb locus to the nearest TSS. Number of bins in each boxplot: n=9962 for B cells p<\/i>\u2009<\u20090.5, n=15877 for B cells p<\/i>\u22650.5, n=12483 for granulocytes p<\/i>\u2009<\u20090.5, n=13345 for granulocytes p<\/i>\u22650.5, n=7337 for NK cells p<\/i>\u2009<\u20090.5, n=18491 for NK cells p<\/i>\u22650.5. Boxplots show 25th percentile, median and 75th percentile, with the whiskers spanning 97% of the data.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 4 Inferring cluster pairs from H3K4me1+H3K27me3 transfers cell type labels.<\/a><\/h3>\n

        (a)<\/b> Histogram of unique fragment cuts per cell. (b)<\/b> Histogram of fraction of unique fragments starting with a “TA” motif. (c)<\/b> UMAP of H3K4me1 sortChIC data, cells colored by cell type. (d)<\/b> Assignment plot showing individual H3K4me1+H3K27me3 cells (represented as dots) assigned to a pair of topics (x-axis labels are H3K4me1 clusters, named by their associated cell type, while y-axis are H3K27me3 clusters). Cells along the diagonal are high-confidence predictions that match a H3K4me1 cluster with a H3K27me3 topics, and are colored by the H3K4me1-derived cell type labels. (e)<\/b> UMAP of H3K4me1+H3K27me3 sortChIC. Cells are colored by their cell type inferred from cluster pairs. Low-confidence predictions are colored in grey. (f, g)<\/b> UMAP representation of H3K4me1 (f) and H3K27me3 (g). Cells are colored by whether the epigenome was generated by single-incubation or by unmixing by scChIX-seq.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 5 Histone modification signal of deconvolved cell types correlates with public H3K4me1 ChIP-seq and H3K27me3 sortChIC ground truth data.<\/a><\/h3>\n

        (a-d)<\/b> Pearson correlation between publicly available H3K4me1 ChIP-seq5<\/a><\/sup> data of purified B cells (a), erythroid (b), granulocytes (c), or NK cells (d) versus H3K4me1 profiles of different cell types derived from scChIX-seq. Single: pseudobulk profiles generated by single incubation, unmixed: pseudobulk profiles deconvolved by scChIX-seq. (e-g)<\/b> Pearson correlation between H3K27me3 sortChIC from FACS-purified B cells (e), granuloytes (f), NK cells (g) versus H3K27me3 sortChIC derived from pseudobulks of whole bone marrow without FACS purification. Single: pseudobulk profiles generated by single incubation, unmixed: pseudobulk profiles deconvolved by scChIX-seq. (h)<\/b> Distribution of assignment probability estimates p<\/i> in the genome for the three cell types. Vertical dotted lines represent cutoffs for p<\/i> to define H3K27me3-specific and H3K4me1-specific regions. p<\/i> is the expected fraction of reads that belong to H3K4me1 in a specific genomic locus. (i)<\/b> Boxplot distributions of GC content for the two classes of regions. (j)<\/b> Boxplot distributions of distance to TSS in the two classes of regions. Distances are measured from the center of the 5 kb locus to the nearest TSS. Boxplots show 25th percentile, median and 75th percentile, with the whiskers spanning 97% of the data.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 6 Re-clustering on B cells reveals heterogeneity within B cells.<\/a><\/h3>\n

        (a)<\/b> UMAP visualization of H3K4me1 and H3K27me3 (single signal and unmixed signal), colored by cell types derived from H3K4me1 and transferred to H3K27me3. Black rectangle indicates the B cell population used to re-cluster in (b,c,d). (b)<\/b> UMAP of pro-B and B cells only. (c,d)<\/b> Projection of H3K4me1 signal of marker genes for pro-B (c) or for differentiated B cells (d). H3K4me1 signal is measured in all cells of the H3K4me1 UMAP (that is both single- and double-incubated have H3K4me1 signal in the H3K4me1 UMAP). Double- (colored) but not single-incubated (grey) cells have H3K4me1 signal in the H3K27me3 UMAP.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 7 H3K4me1 and H3K27me3 signal during neutrophil maturation.<\/a><\/h3>\n

        (a)<\/b> UMAP visualization of H3K4me1 and H3K27me3, lines join H3K4me1 and H3K27me3 UMAPs of double-incubated neutrophils. Heterogeneity within neutrophils are colored as neutrophil pseudotime. (b)<\/b> H3K4me1 and H3K27me3 modification levels at the Retnlg<\/i> (a mature neutrophil marker gene) locus along neutrophil pseudotime. (c)<\/b> H3K4me1 and H3K27me3 modification levels at the Hoxa<\/i> along neutrophil pseudotime. (d)<\/b> UMAP of H3K27me3 signal across single cells colored by weights of a topic containing high H3K27me3 levels at many Hox<\/i> and developmental gene loci (Hox<\/i> topic). (e)<\/b> Topic weights of the top 150 genes associated with loci in the Hox<\/i> topic for H3K27me3. (f)<\/b> Neutrophil mRNA abundance of genes in the Hox<\/i> topic compared to other genes derived from publicly available scRNA-seq data25<\/a><\/sup>. Number of genes per boxplot: n=17986 for All Genes, n=127 for genes in the Hox topic. Boxplots show 25th percentile, median and 75th percentile, with the whiskers spanning 97% of the data.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 8 Cell typing mouse organogenesis dataset using H3K36me3 using marker genes.<\/a><\/h3>\n

        (a)<\/b> Histogram of unique fragment cuts per cell. (b)<\/b> Histogram of fraction of unique fragments starting with a “TA” motif. (c-l)<\/b> Genome browser plots of cell type-specific H3K36me3 loci showing pseudobulk CPM signals (colored lines, top) and cut locations of individual cells (bottom, black marks). Cells are ordered by cell type (color-coded on the left).<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 9 H3K9me3-specific regions across cell types.<\/a><\/h3>\n

        (a)<\/b> Heatmap of H3K36me3 signal for the top 250 H3K36me3-specific loci (rows) across cell types (columns). (b)<\/b> Heatmap of mRNA abundances for the genes associated with the H3K36me3-specific loci in (a) across pseudobulks. Data processed from publicly available scRNA-seq data from Cao et al.42<\/a><\/sup>. (c)<\/b> Heatmap of H3K9me3 signal for the same top 250 H3K36me3-specific loci as in (a). The H3K36me3 and H3K9me3 heatmaps are mean-centered and scaled using a common mean and standard deviation calculated across both marks. (d)<\/b> Heatmap of H3K9me3 signal across pseudobulks at H3K9me3-variable loci. (e)<\/b> Relative mRNA abundances42<\/a><\/sup> at n=364 genes associated with erythroblast-repressed loci across nine cell types. (f)<\/b> mRNA abundance of an erythroblast-repressed gene, Nell2<\/i>, across pseudobulks. (g)<\/b> Genome browser plot of around the Nell2<\/i> locus, an erythroblast-specific region for H3K9me3. Top of plot is pseudobulk H3K9me3 CPM signals, below are cut locations of individual cells (black marks). Cells are ordered by cell type (color-coded as in heatmaps). (h, i)<\/b> Total unique fragments across cell types for single-incubated cells for H3K36me3 (h) and H3K9me3 (i), showing that the variability of the number of cuts across cells can span orders of magnitude. Number of single-incubated H3K36me3 cells for each boxplot: n=154 erythroid, n=36 white blood cells, n=60 endothelial, n=250 neural tube progenitors, n=272 neurons, n=58 Schwann cell precursors, n=154 epithelial, n=570 mesenchymal progenitors, n=160 cardiomyocytes. For H3K9me3: n=207 erythroid, n=26 white blood cells, n=736 non-blood cell types. Boxplots in (e), (h), (i) show 25th percentile, median and 75th percentile, with the whiskers spanning 97% of the data.<\/p>\n<\/div>\n

        \n

        Extended Data Fig. 10 Distinct dynamics of H3K4me1 and H3K36me3 during macrophage in vitro differentiation.<\/a><\/h3>\n

        (a)<\/b> Density plots of total number of cuts across cells for H3K4me1, H3K36me3, and H3K4me1+H3K36me3 labeled cells. (b)<\/b> Density plots of fraction of cuts starting with a TA motif across cells for H3K4me1, H3K36me3, and H3K4me1+H3K36me3 labeled cells. (c)<\/b> Genome-browser plot around gene body of Mertk<\/i>, a macrophage-specific gene. Tracks are bigwigs from pseudobulks averaged across the time course. (d)<\/b> Log2 fold change estimates along pseudotime on gene bodies in the genome. Colored dots are considered significant (log2 fold change in H3K36me3 > 3.5, zscore in both H3K36me3 and H3K4me1 > 2) and used for chromatin velocity estimates. (e, f)<\/b> UMAP of H3K4me1 (e) and H3K36me3 (f) of single-incubated and deconvolved cells showing intermingling of the two types of cells. (g)<\/b> Examples of H3K4me1 and H3K36me3 for an upregulated (above) and downregulated (below) gene along pseudotime. (h)<\/b> Histogram of estimates of the rate constant \u03b3<\/i> for the 209 dynamic genes highlighted in (d).<\/p>\n<\/div>\n<\/div>\n<\/div>\n

        \n

        Supplementary information<\/h2>\n<\/div>\n