
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae426
bbae426
Problem Solving Protocol
AcademicSubjects/SCI01060
SnapHiC-G: identifying long-range enhancer–promoter interactions from single-cell Hi-C data via a global background model
https://orcid.org/0000-0003-2278-8946
Liu Weifang Department of Biostatistics, University of North Carolina at Chapel Hill, 135 Dauer Drive, Chapel Hill, NC 27599, United States

Zhong Wujuan Biostatistics and Research Decision Sciences, Merck & Co., Inc., 126 East Lincoln Ave, Rahway, New Jersey 07065, United States

Giusti-Rodríguez Paola Department of Psychiatry, University of Florida, 1149 Newel Dr., Gainesville, FL 32611, United States

Jiang Zhiyun Department of Genetics, University of North Carolina at Chapel Hill, 120 Mason Farm Road, Chapel Hill, NC 27599, United States

Wang Geoffery W Department of Biostatistics, University of North Carolina at Chapel Hill, 135 Dauer Drive, Chapel Hill, NC 27599, United States

Sun Huaigu Department of Genetics, University of North Carolina at Chapel Hill, 120 Mason Farm Road, Chapel Hill, NC 27599, United States

https://orcid.org/0000-0003-0987-2916
Hu Ming Department of Quantitative Health Sciences, Lerner Research Institute, Cleveland Clinic Foundation, 9500 Euclid Avenue, Cleveland, OH 44196, United States

https://orcid.org/0000-0002-9275-4189
Li Yun Department of Biostatistics, University of North Carolina at Chapel Hill, 135 Dauer Drive, Chapel Hill, NC 27599, United States
Department of Genetics, University of North Carolina at Chapel Hill, 120 Mason Farm Road, Chapel Hill, NC 27599, United States
Department of Computer Science, University of North Carolina at Chapel Hill, 201 S. Columbia St, Chapel Hill, NC 27599, United States

Corresponding authors. Department of Quantitative Health Sciences, Lerner Research Institute, Cleveland Clinic Foundation, Cleveland, OH, USA. E-mail: hum@ccf.org; Department of Biostatistics, Department of Genetics, Department of Computer Science, University North Carolina, Chapel Hill, NC 27599, USA. E-mail: yunli@med.unc.edu
Weifang Liu and Wujuan Zhong contributed equally to this work.

9 2024
02 9 2024
02 9 2024
25 5 bbae42620 1 2024
05 7 2024
13 8 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Harnessing the power of single-cell genomics technologies, single-cell Hi-C (scHi-C) and its derived technologies provide powerful tools to measure spatial proximity between regulatory elements and their target genes in individual cells. Using a global background model, we propose SnapHiC-G, a computational method, to identify long-range enhancer–promoter interactions from scHi-C data. We applied SnapHiC-G to scHi-C datasets generated from mouse embryonic stem cells and human brain cortical cells. SnapHiC-G achieved high sensitivity in identifying long-range enhancer–promoter interactions. Moreover, SnapHiC-G can identify putative target genes for noncoding genome-wide association study (GWAS) variants, and the genetic heritability of neuropsychiatric diseases is enriched for single-nucleotide polymorphisms (SNPs) within SnapHiC-G-identified interactions in a cell-type-specific manner. In sum, SnapHiC-G is a powerful tool for characterizing cell-type-specific enhancer–promoter interactions from complex tissues and can facilitate the discovery of chromatin interactions important for gene regulation in biologically relevant cell types.

Hi-C
enhancer–promoter interactions
single cell
epigenetics
GWAS
National Institutes of Health 10.13039/100000002 R35HG011922 U01DA052713 R01MH125236 P50HD103573 R01NR019245
==== Body
pmcIntroduction

Chromatin spatial organization plays an essential role in genome function and gene regulation [1–3]. Genome-wide chromatin conformation capture technologies like Hi-C have improved our understanding of the 3D genome structure and aided the discovery of chromatin features at various scales, such as chromosome territories, A/B compartments, topologically associating domains (TADs), and chromatin loops. In bulk Hi-C analysis, millions of cells are pooled together to generate a population-averaged chromatin contact map. However, chromatin features obtained in a large population of cells are not representative of individual cells from a complex mixture [4, 5], since the dynamics of chromatin conformation in individual cells are masked. Although larger-scale organizational features of the genome are evolutionarily conserved and highly reproducible at the population level [6], direct evidence from single-cell imaging and sequencing-based technologies reveals considerable cell-to-cell variability in chromatin folding features at smaller units such as chromatin loops and enhancer–promoter interactions [7–13], suggesting that genome organization is remarkably more variable than expected [14]. In particular, chromatin interactions identified in bulk Hi-C data are typically present in only a relatively small proportion (10%–30%) of the cell population at any given time [8, 9], indicating a high degree of cell-level heterogeneity, consistent with the stochastic gene expression widely observed in mammals [14]. Therefore, studying chromatin features in single cells can provide new insights into genome function and gene regulation.

Recent advances in single-cell genomics enable profiling of 3D genome structures in single cells at unprecedented throughput and resolution [7, 11, 12, 15–21]. Single-cell Hi-C (scHi-C) adopts the bulk Hi-C protocol with an extra step to isolate or barcode single cells [15]. Like bulk Hi-C, the standard scHi-C data processing workflow includes read alignment, deduplication, filtering of contacts and cells, and binning to generate the final contact matrix ready for data analysis [15, 16]. Specifically, chromatin contacts for every cell are represented by a symmetric matrix with entries representing contact frequencies between fixed-length genomic regions, referred to as bin pairs. Analysis of scHi-C data includes 3D architecture reconstruction and feature extraction (enhancer–promoter interactions, chromatin loops, TADs, A/B compartments) for single cells [15, 16]. Due to the extreme sparsity in scHi-C data (over 99% zeros in contact matrices at kilobase resolution) and limited capture efficiency, imputation of missing contacts is essential to enhance the signal-to-noise ratio and subsequently allow for the detection of chromatin features. The first scHi-C data imputation method was proposed by Zhou et al. [22], which applies the random walk with restart (RWR) algorithm. An alternative method, Higashi [23, 24], models both cells and bins as nodes and then uses a hypergraph neural network to learn edges between these nodes, considering the similarities and variability between cells instead of imputing each cell separately. More recently, scVI-3D [25] was developed to model interaction frequencies by a deep degenerative model with a zero-inflation component to account for sparsity in scHi-C data. scVI-3D can impute sparse scHi-C contacts and account for the impact of genomic distance, library size, and batch effects. Furthermore, with imputed data, latent representations of scHi-C can be learned to reconstruct 3D genome structures like TADs [26].

While computational methods tailored for scHi-C data are still under development, many approaches have been proposed to detect chromatin loops or statistically significant long-range chromatin interactions from bulk Hi-C data. They can be broadly grouped into two classes, namely, global background–based methods and local background–based methods [27–31]. Global background–based methods fit a global statistical model based on 1D genomic distance and assign P-values to each bin pair in the contact matrix by comparing the observed contact frequency to the expected contact frequency under the global background. In contrast, local background–based methods identify peaks in the Hi-C contact map that are local maxima with respect to their neighboring bin pairs. Our recent work SnapHiC and SnapHiC2 [32, 33] perform data imputation by RWR first and then combine both global and local background models to identify chromatin loops from single cells of the same cell type. While being the first and only existing method to detect chromatin loops from scHi-C data, most SnapHiC-identified chromatin loops are CCCTC-binding factor (CTCF)-anchored structural loops due to the use of the local background model, while the sensitivity to identify enhancer–promoter interactions is relatively low. To the best of our knowledge, no method exists for explicitly identifying enhancer–promoter interactions from scHi-C data.

To fill this gap, we propose SnapHiC-G, a new computational approach based solely on a global background model to identify long-range enhancer–promoter interactions from scHi-C data. We applied SnapHiC-G to re-analyze scHi-C datasets generated from mouse embryonic stem cells (mESCs) and human brain cortical cells. We showed that SnapHiC-G outperformed SnapHiC and existing methods designed for bulk Hi-C data, achieving higher sensitivity with comparable precision in identifying long-range enhancer–promoter interactions.

Results

Overview of the SnapHiC-G algorithm

SnapHiC-G identifies enhancer–promoter interactions based on both scHi-C data and epigenetic annotations. The algorithm consists of four components: (i) imputing chromatin contact probabilities in each single cell, (ii) distance-stratified normalization of imputed contact probabilities, (iii) filtering candidate bin pairs, and (iv) identifying statistically significant long-range enhancer–promoter interactions (Fig. S1). Given single-cell contact matrices, SnapHiC-G first applies the RWR algorithm to generate imputed contact probability matrices for each single cell using a sliding window approach [33] (SnapHiC-G algorithm). Next, the imputed contact probabilities are converted into distance-stratified Z-scores to account for the dependence between contact probability and the 1D genomic distance between two bins. To identify enhancer–promoter interactions, SnapHiC-G filters bin pairs based on transcript start sites (TSSs) and available epigenetic annotations (e.g. H3K4me3 and H3K27ac ChIP-seq peaks or ATAC-seq peaks) to obtain a set of candidate bin pairs that span between gene promoters and cis-regulatory elements (SnapHiC-G algorithm). SnapHiC-G then defines enhancer–promoter interactions based on the global background by applying a one-sample t-test for each tested bin pair across all single cells belonging to the same cell type. Specifically, for each bin pair, a one-sided hypothesis test is conducted where the null hypothesis states that the bin pair’s average normalized contact probability across all single cells equals zero. The alternative hypothesis states that the average normalized contact probability is greater than zero. By default, bin pairs with false discovery rate (FDR)<0.1 and t-statistics>3 are identified as significant enhancer–promoter interactions (SnapHiC-G algorithm). In our analysis, all scHi-C data are binned into 10Kb resolution unless stated otherwise.

Benchmarking with mouse embryonic stem cells

We applied SnapHiC-G, SnapHiC [32], Chromosight [27], FitHiC2 [28], FastHiC [29], HiC-ACT [30], and HiC-DC+ [31] on single-cell Hi-C data generated from 742 mESCs [18], where the latter five are methods designed for bulk Hi-C data. We aggregated single-cell Hi-C data for all cells as a pseudo-bulk Hi-C sample as input for bulk Hi-C methods (Identification of loops/interactions using  other Hi-C methods). To benchmark against other methods, a reference list of significant interactions was constructed using HiCCUPS-identified loops [34] from deeply sequenced bulk Hi-C data [35] and model-based analysis of PLAC-seq (MAPS) and HiChIP-identified significant interactions from H3K4me3 PLAC-seq [36], cohesin HiChIP [37], and H3K27ac HiChIP data [38]. We took the union of all reference interactions and kept only bin pairs with a genomic distance between 20Kb and 1Mb for evaluation. Since interactions called by SnapHiC-G were filtered with external epigenetic data, to ensure a fair comparison and demonstrate its potential to detect enhancer–promoter interactions, we applied the same filtering steps to results from other methods.

When applying to the complete set of 742 mESCs, SnapHiC-G identified notably more significant enhancer–promoter interactions than other methods and reached a genome-wide power of 80%, which means that among all the 38 588 interactions in the reference list, SnapHiC-G-identified interactions recovered 80% of them (Fig. 1A; Table 1). Due to extreme data sparsity, bulk Hi-C methods missed most reference interactions without imputing single-cell contact probabilities. FitHiC2 performed the best among other methods, calling 3476 interactions with a genome-wide power of 16%. FastHiC identified more interactions than FitHiC2 with a lower genome-wide power of 12%. As expected, SnapHiC identified few enhancer–promoter interactions due to the use of the local background model. Although already tailored for sparse scHi-C data, SnapHiC performed similarly to HiC-ACT, a global background method. HiC-DC+ identified the fewest interactions among all methods with the lowest genome-wide power. Moreover, we benchmarked SnapHiC-G identified interactions with clustered regularly interspaced short palindromic repeats (CRISPR)-validated enhancer–promoter interactions in mESCs. Results from 742 mESCs showed that SnapHiC-G successfully detected previously verified long-range enhancer–promoter interactions at Sox2, Atpif1, Phactr4, and Med13l loci [39–42], while SnapHiC and FitHiC2 were only able to detect the interaction at Sox2 locus.

Figure 1 Power curves with (A) 742 mESCs and (B) 100 mESCs. Interactions were ranked by significance on the x-axis and power was evaluated with the corresponding number of top interactions. The sub-figures at the lower right-hand corner are zoomed-in views of the top 10 000 interactions for each method.

Table 1 Genome-wide power and number of interactions called for mESCs and three brain cell types. Each row represents a method, and each column represents a data set. # interactions: number of significant interactions called by each method. Power: the proportion of significant interactions that overlapped with the reference list. Results are shown for all genome-wide significant interactions. The number before each cell type represents the number of cells in each data set.

	742 mESCs	100 mESCs	261 L2/3 neurons	323 microglia	1038 oligodendrocytes	
Method	# interactions	Power	# interactions	Power	# interactions	Power	# interactions	Power	# interactions	Power	
HiC-DC+	582	0.050	59	0.007	4359	0.283	4321	0.245	9881	0.466	
FastHiC	3988	0.116	3576	0.073	6585	0.207	7852	0.245	15 583	0.492	
HiC-ACT	1365	0.074	53	0.005	1471	0.159	1881	0.172	9166	0.575	
FitHiC2	3476	0.160	67	0.005	4019	0.298	4834	0.310	18 708	0.723	
Chromosight	1980	0.061	1935	0.055	3068	0.108	2789	0.081	3554	0.118	
SnapHiC	1070	0.076	265	0.028	2918	0.203	1766	0.147	3605	0.242	
SnapHiC-G	121 469	0.793	81 830	0.607	236 802	0.879	165 898	0.781	653 890	0.873	

Since the number of input cells is critical in scHi-C data analysis, we assessed whether SnapHiC-G could retain its performance with fewer cells. Among all 742 mESCs, 100 mESCs were randomly selected, and the same performance evaluation of the seven methods mentioned above was repeated. As shown in Fig. 1B and Table 1, all methods had reduced power with 100 mESCs, while SnapHiC-G was least affected by the number of cells and showed more significant power gain over other methods compared with results from 742 cells. With only 100 mESCs, SnapHiC-G retained a genome-wide power of 61%, while all the other methods had genome-wide power below 10%. FastHiC still performed the best among others, with 3576 significant interactions identified, comparable to the complete data, but the genome-wide power reduced almost by half to 7%. On the other hand, SnapHiC reached a 3% genome-wide power with 265 loops called; however, it is still more sensitive than other bulk Hi-C methods.

Due to the relatively large number of significant interactions detected by SnapHiC-G, we evaluated the precision of identified interactions across the seven methods with the same reference list. To control for the number of bin pairs compared, we ranked identified interactions from each method by their significance (i.e. P-value for HiC-ACT or FDR for FitHiC2, HiC-DC+, SnapHiC), posterior probability (FastHiC), or correlation score (Chromosight) and calculated precision for the top 1000, 2000, 5000, and 10 000 interactions. SnapHiC-G showed a comparable or better performance in terms of precision among the most significant interactions, even with a much larger number of interactions called (Fig. 2). For example, with 742 mESCs, SnapHiC-G attained a precision of 0.93 for the top 1000 interactions, which was comparable with HiC-DC+ (0.96) and FastHiC (0.96) and substantially higher than HiC-ACT (0.73), FitHiC2 (0.73), SnapHiC (0.72), and Chromosight (0.46). With 100 mESCs, SnapHiC-G had a precision of 0.60–0.76 for the top 1000 to 10 000 interactions. We acknowledge that the reference list is not the gold standard and may still miss many true enhancer–promoter interactions; therefore, it is particularly challenging to define true negatives. As an alternative, we benchmarked the specificity of SnapHiC-G-identified interactions using a number of datasets reporting experimentally tested enhancer–promoter connections in mESCs collected by Fulco et al. [43]. Specifically, we examined enhancer–promoter pairs that did not exhibit a significant change in gene expression after CRISPR perturbation experiments and defined these as true negatives. Out of the 200 unique true negative enhancer–promoter pairs, only 37 were identified as significant by SnapHiC-G, resulting in a specificity of 81.5%. This suggests that SnapHiC-G demonstrates reasonable performance in terms of specificity. We performed similar evaluation with three cell types from human brain cortical cells (oligodendrocytes, microglia, and L2/3 neurons) and observed that SnapHiC-G has higher sensitivity compared with other methods (Supplementary Text; Figs S2 and S3).

Figure 2 Precision bar plots with (A) 742 mESCs and (B) 100 mESCs. Results shown are for the top 1000, 2000, 5000, and 10 000 interactions ranked by significance. Some bars are missing because the number exceeds the number of interactions called by that method.

Although SnapHiC-G is designed to retrieve enhancer–promoter interactions, the identified interactions contained a small proportion of structural interactions. Using 742 mESCs, we compared the fraction of structural interactions identified by SnapHiC-G, SnapHiC, and FitHiC2, both by overlapping and exclusive interactions called by each method. A structural interaction is defined as those where both anchor bins overlap with a CTCF ChIP-seq peak region [44]. Results show that SnapHiC-G-specific interactions and FitHiC2-specific interactions had a similar proportion (~21%) of structural interactions, which are lower than interactions identified by two or three methods (28%–50%) (Fig. S4; Table S1).

In addition to the analysis of 10Kb resolution data, the application of SnapHiC-G to 5Kb resolution data for 100 mESCs yielded compelling results, showcasing a precision of 0.82 for the top 1000 interactions. This underscores SnapHiC-G’s capability in analyzing 5Kb-resolution scHi-C data with high precision and reasonable detection power (Table S2). Taken together, we have shown that SnapHiC-G achieved much higher sensitivity than other methods while maintaining precision among top enhancer–promoter interactions even with a small number of cells.

Enrichment of brain expression quantitative trait locus (eQTL)-TSS pairs in human cortical cells

To evaluate the functional characteristics of SnapHiC-G identified enhancer–promoter interactions, we re-analyzed single-nucleus methyl-3C-seq (sn-m3c-seq) data from 2869 human prefrontal cortical cells [45]. For sn-m3c-seq data, cell types were classified based on DNA methylome as described in the original study [45]. We applied SnapHiC-G to four major cell types [astrocytes (n = 338), microglia (n = 323), oligodendrocytes (n = 1038), and L2/3 neurons (n = 261)] separately, where single cells from the same cell type were pooled (Fig. 3A; Table 1; Table S3).

Figure 3 (A) Enrichment of SnapHiC-G-identified interactions for each cell type in “brain cell type”, CommonMind, and GTEx liver eQTL-TSS pairs. Here “brain cell type” refers to the eQTL data from Bryois et al., the squares denote point estimates for odds ratios (ORs) and the error bars denote 95% confidence intervals for ORs. OR was defined as the ratio of the odds of SnapHiC-G-identified interactions overlapping with eQTL-TSS bin pairs over the odds of pseudo bin pairs overlapping with eQTL-TSS bin pairs. SnapHiC-G interactions were identified from four major cell types (astrocytes (Astro) [n = 338 cells], microglia (MG) [n = 323], oligodendrocytes (ODC) [n = 1038], and L2/3 neurons (L2/3) [n = 261]). (B) UpSet plot for SnapHiC-G-identified interactions. Y-axis is the number of exact overlapped interactions between four cell types (astrocytes (Astro), microglia (MG), oligodendrocytes (ODC), and L2/3 neurons (L2/3) (Definition of exact overlapped enhancer–promoter interactions). The interactions were identified from 261 cells from each cell type. The UpSet figure was generated using the R package UpSetR (Web Resources for figure generation).

We created control bin pairs to compute the enrichment of SnapHiC-G-identified enhancer–promoter interactions as eQTL-TSS of eGene pairs for each cell type. Specifically, we generated a pseudo bin pair for each significant enhancer–promoter interaction by retaining the bin with the promoter as the center but flipping the other bin to the opposite side of the center (Enrichment of eQTL-TSS pairs). Additionally, we collected eQTL data from three sources, including “brain cell type” eQTL data from Bryois et al. [46], the CommonMind Consortium (CMC) brain eQTL data [47], and (as a control sample) GTEx consortium v7 liver eQTL data [48]. Here the “brain cell type” eQTL data were for separate brain cell types released from the Netherlands Brain Bank (NBB), the MS UK Tissue Bank (UKTB), and the Edinburgh Brain Bank (EBB) [46].

We found that the odds of a bin pair to be an eQTL-TSS pair were significantly higher for significant interactions than pseudo bin pairs by overlapping both with true eQTL-TSS bin pairs (Fig. 3A). The wide confidence intervals (CI) of the odds ratios (ORs) calculated from the Bryois et al. data [46] were due to a smaller sample size to detect eQTLs compared with the CMC brain eQTL data (192 versus 467 brain samples, respectively). As expected, SnapHiC-G-identified interactions were more enriched in brain eQTLs than in liver eQTLs for most cell types, and interactions were more enriched in eQTLs identified from matched brain cell types than from un-matched brain cell types. For example, the ORs of eQTLs from Bryois et al., CMC brain, and GTEx liver in microglia were 1.41 (95% CI: 1.07–1.87), 1.22 (95% CI: 1.18–1.26), and 1.05 (95% CI: 0.93–1.18), respectively. To ensure a comparable number of interactions were being compared, we selected top interactions from SnapHiC-G to match the number of loops identified by SnapHiC. We then repeated the enrichment analysis and observed consistent results (Fig. S5). These results demonstrate the functional importance of SnapHiC-G identified enhancer–promoter interactions in relevant cell types and tissues. We also compared SnapHiC-G with SnapHiC and FastHiC regarding the enrichment of eQTL-TSS pairs. Our findings revealed that interactions identified by SnapHiC-G exhibit much narrower 95% CI for ORs compared to interactions identified by SnapHiC and FastHiC (Fig. S6).

Gene expression patterns of cell-type-specific enhancer–promoter interactions

To further assess the biological relevance of SnapHiC-G-identified enhancer–promoter interactions, we evaluated gene expression patterns in cell-type-specific interactions from the four brain cell types. To adjust for the impact of different numbers of cells across different cell types, we downsampled astrocytes, microglia, and oligodendrocytes to 261 cells each to match the number of L2/3 neuron cells. We then identified cell-type-specific interactions with SnapHiC-G from the downsampled data (Down-sampling of scHi-C data). We defined cell-type-specific interactions as those present exclusively in one cell type, not including those shared by two or more cell types (Fig. 3B; Table S4; Definition of overlapped enhancer–promoter interactions). Additionally, we assessed gene expression patterns associated with SnapHiC-identified chromatin loops to facilitate a comprehensive comparison with SnapHiC-G. Notably, the same set of 261 cells was used in both SnapHiC-G and SnapHiC analyses for consistency. Next, we selected genes with promoters overlapping with cell-type-specific interactions and extracted gene expression levels from the corresponding cell types [49]. As shown in Fig. 4A and B, gene expression levels were higher for genes from cell-type-specific enhancer–promoter interactions in the corresponding cell type, for both SnapHiC-G and SnapHiC. In most scenarios, SnapHiC-G exhibits more significant P-values than SnapHiC when testing the differences in gene expression levels between cell types. This observation can be attributed to the larger number of interactions identified by SnapHiC-G as compared to SnapHiC. Moreover, after removing genes that overlapped with cell-type-specific enhancer–promoter interactions detected by SnapHiC, most (>83%) of the SnapHiC-G identified genes remained, highlighting the substantial additional information from SnapHiC-G. We also observed increased expression of these genes in a cell-type-specific manner (Fig. 4C). Similar to the eQTL-TSS enrichment analysis, we selected top interactions from SnapHiC-G to match the number of loops identified by SnapHiC and repeated the cell-type-specific gene expression analysis and observed similar patterns (Fig. S7). Furthermore, we evaluated the enrichment of cell-type-specific gene expression in SnapHiC-G interactions and observed that astrocyte-specific SnapHiC-G interactions showed higher enrichment for genes specifically expressed in astrocytes than the genes specifically expressed in microglia, oligodendrocytes, or neurons (Fig. S8). We had similar observations for the enrichment of microglia-, oligodendrocytes-, and L2/3 neurons-specific SnapHiC-G interactions. These results add another line of evidence that predicted enhancer–promoter interactions could provide valuable information in a cell-type-specific manner, complementary to the eQTL enrichment analysis.

Figure 4 Cell-type-specific gene expression: Violin plots of RNA-seq based expression levels (log2(FPKM+1)) for selected genes. (A) Genes overlapping with cell-type-specific SnapHiC-G enhancer–promoter interactions. (B) Genes overlapping with cell-type-specific SnapHiC enhancer–promoter interactions. (C) Genes overlapping with cell-type-specific SnapHiC-G enhancer–promoter interactions, where those overlapping with cell-type-specific enhancer–promoter interactions detected by SnapHiC were subsequently removed. P-values were calculated from paired Wilcoxon signed-rank tests. Gene expression outliers for each cell type were removed for visualization. For both SnapHiC-G and SnapHiC, the interactions were identified from the same 261 cells from each cell type.

Assigning GWAS variants to putative target genes

With over 90% of GWAS variants associated with human complex diseases and traits residing in noncoding regions yet enriched in cis-regulatory elements (e.g. promoters, enhancers, silencers, and insulators), enhancer–promoter interactions have the potential to prioritize disease-relevant genes for noncoding variants, particularly those in close spatial proximity that are far away in the 1D genomic distance with the promoter of their target genes. To assign putative target genes to noncoding GWAS variants based on predicted enhancer–promoter interactions in brain cell types, we collected the latest GWAS summary statistics for eight neurodevelopmental and neurodegenerative disorders: Alzheimer’s disease (AD) [50], attention-deficit hyperactivity disorder (ADHD) [51], autism spectrum disorders (ASDs) [52], bipolar disorder (BP) [53], schizophrenia (SCZ) [54], Parkinson’s disease (PD) [55], major depressive disorder (MDD) [56], and neuroticism (NEU) [57] and two complex traits: educational attainment (EDU) [58] and intelligence quotient (IQ) [59]. We again focused on cell-type-specific enhancer–promoter interactions to predict target genes in a cell-type-specific manner using downsampled sn-m3C-seq data. Specifically, SnapHiC-G identified 137 418, 154 261, 157 181, and 236 802 enhancer–promoter interactions in astrocytes, microglia, oligodendrocytes, and L2/3 neurons, respectively (Table 2). After excluding interactions shared among cell types, 14 440, 39 220, 23 853, and 65 819 cell-type-specific interactions were left correspondingly. To facilitate the interpretation of the GWAS variants, we focused on noncoding GWAS variants that reside in an active enhancer region in astrocytes, microglia, oligodendrocytes, or L2/3 neurons [60]. When matching GWAS variants and SnapHiC-G results, for each cell-type-specific enhancer–promoter interaction, we required that one bin contains GWAS variant(s) and the other bin overlaps with a gene’s TSS, and we annotated this gene as the putative target gene. Furthermore, we required that the corresponding gene is highly expressed [fragments per kilobase of transcript per million mapped reads (FPKM)>1] in this cell type and lowly expressed (FPKM ≤1) in the other three cell types. We found 35, 82, 7, and 98 matched enhancer–promoter interactions (222 in total) for astrocytes, microglia, oligodendrocytes, and L2/3 neurons, respectively, and resolved over 600 SNP–disease associations (Table 2). Moreover, the average number of target genes for each variant was close to 1, much smaller than the number of nearby genes (+/−1Mb region), ranging from 25 to 83. For example, in astrocytes, the average number of target genes and nearby genes per variant were 1.1 and 38.3, respectively. To further evaluate the functional importance of matched enhancer–promoter interactions, we collected fine-mapped GWAS SNPs for AD from Schwartzentruber et al. [50, 61] and for SCZ from Trubetskoy et al. [62]. We identified cell-type-specific interactions from 261 L2/3 neurons and 323 microglia using SnapHiC-G, SnapHiC, and FitHiC2 and asked how many of the matched GWAS variants from each method are fine-mapped. Out of the 331 unique matched AD GWAS variants identified by SnapHiC-G in microglia, 29 were fine-mapped; out of the 18 unique matched SCZ GWAS variants identified by SnapHiC-G in L2/3 neurons, 9 were fine-mapped. In comparison, SnapHiC did not map any fine-mapped SNPs to target genes in microglia or L2/3 neurons, and FitHiC2 was only able to map one fine-mapped variant for SCZ in L2/3 neurons. These results showed that the reported interactions are indeed functionally relevant and target genes of noncoding variants can be pinpointed in a cell-type-specific manner by integrating SnapHiC-G identified enhancer–promoter interactions with GWAS results.

Table 2 Summary of SnapHiC-G-identified enhancer–promoter interactions for down-sampled 261 cells from each brain cell type. SnapHiC-G interactions: number of enhancer–promoter interactions identified from SnapHiC-G; cell-type-specific interactions: number of enhancer–promoter interactions identified only in this cell type and not in any other three cell types; matched interactions: number of cell-type-specific enhancer–promoter interactions with one bin containing the GWAS SNP residing in an active enhancer while the other bin overlapping with a gene’s TSS and the corresponding gene is highly expressed (FPKM>1) in this cell type and lowly expressed (FPKM ≤1) in other three cell types; unique GWAS SNPs: number of unique SNPs contained in the matched enhancer–promoter interactions; SNP disease associations: number of GWAS disease−SNP associations in the matched enhancer–promoter interactions; average number of target genes per GWAS SNP: average number of targeted genes for each GWAS SNP based on the matched enhancer–promoter interactions; average number of (+/− 1Mb) genes per GWAS SNP: average number of nearby genes (within 1Mb) for each GWAS SNP. The interactions were identified from 261 cells from each cell type.

Cell type	SnapHiC-G interactions	Cell-type-specific interactions	Matched interactions	Unique GWAS SNPs	SNP–disease Associations	Average number of target genes per GWAS SNP	Average number of (+/− 1Mb) genes per GWAS SNP	
Astrocytes	137 418	14 440	35	72	72	1.07	38.32	
Microglia	154 261	39 220	82	279	288	1.03	82.86	
Oligodendrocytes	157 181	23 853	7	22	39	1.00	46.91	
L2/3 neurons	236 802	65 819	98	181	202	1.02	24.61	

Examples of cell-type-specific enhancer–promoter interactions

From the 222 matched cell-type-specific enhancer–promoter interactions mentioned in the previous section, we were able to map a GWAS variant residing in an active enhancer with the promoter of a cell-type-specific gene (Table 2). Notably, most of these interactions were not identified by SnapHiC with the local background approach, as shown in Fig. 4B. We demonstrate how SnapHiC-G-identified enhancer–promoter interactions can elucidate functional genes of GWAS loci in relevant cell types with a few examples.

The first example locates at a locus on chromosome 8, showing cell-type-specific interactions in L2/3 neurons and astrocytes (Fig. 5; Fig. S9). Specifically, an SCZ-associated GWAS SNP rs2565064 (chr8: 27 327 841) interacts with the promoter of PNOC in neurons, while another SCZ-associated SNP rs28541694 (chr8: 27 462 008) interacts with the promoter of ZNF395 in astrocytes. In addition, 10 AD-associated SNPs located in the chr8:27.4Mb–27.5Mb region are also connected to the promoter region of ZNF395 in astrocytes, while none of the AD-associated SNPs interact with the PNOC promoter. These two genes also show consistent cell-type-specific gene expression patterns, with PNOC being highly expressed in neurons (FPKM = 4.75 in neurons versus ≤1 in the other three cell types) and ZNF395 being highly expressed in astrocytes (FPKM = 2.26 in astrocytes versus ≤1 in the other three cell types). Moreover, rs2565064 resides in a neuron-specific enhancer, and SNPs interacting with ZNF395 reside in an astrocyte-specific enhancer. These results indicated that PNOC is the putative target gene for SCZ-associated SNP rs2565064 in neurons and ZNF395 is the putative target gene for both SCZ- and AD-associated SNPs in this 27.4Mb–27.5Mb region on chromosome 8. PNOC is primarily transcribed in the brain and spinal cord in the central nervous system [63] and encodes the precursor for bioactive neuropeptides that influence a broad range of physiological roles, including memory, learning, and neuronal development, fear, anxiety, and sleep [64]. PNOC was also associated with post-traumatic stress disorder (PTSD), whose transcriptome significantly correlated with SCZ [65, 66]. On the other hand, ZNF395, a gene involved in inflammation and cancer progression [67], is upregulated in the SCZ network from the integrative network analysis [68].

Figure 5 An illustrative example at the PNOC-ZNF395 locus on chromosome 8 with neuron- and astrocyte-specific interactions. The top panel shows the gene track. The middle panels show RNA-seq, H3K27ac, and ATAC-seq tracks for the four brain cell types. The bottom panels show SCZ and AD GWAS SNPs, enhancer regions in astrocytes and neurons, and cell-type-specific interactions identified by SnapHiC-G but not SnapHiC. Astrocyte-specific interactions link the promoter region of ZNF395 (highlighted bar on the right) to a SCZ-associated SNP rs28541694 and 10 AD-associated SNPs located in the chr8:27.4Mb–27.5Mb locus (highlighted bar on the left); a neuron-specific interaction links the promoter region of PNOC (highlighted bar on the right) to another SCZ-association SNP rs2565064 (highlighted bar on the left). Both genes highlighted showed celltype-specific gene expression in corresponding cell types. The anchors of these interactions also showed stronger H3K27ac ChIP-Seq and ATAC-seq signals in matched cell types.

Next, we focused on a microglia-specific interaction between the promoter of ARPC1B, which is highly expressed in microglia (FPKM = 7.12 in microglia versus ≤ 1 in the other three cell types), and an AD-GWAS locus located in an active enhancer on chromosome 7 at 99.7Mb (Fig. S10). Our results showed that ARPC1B is the predicted target gene for this AD-associated GWAS variant rs1880949, consistent with prior findings that ARPC1B was active in microglia in AD patients but not in healthy controls [69]. At the same locus, SnapHiC-G detected another microglia-specific enhancer–promoter interaction between an active enhancer containing an EDU-associated GWAS variant rs10241492 (chr7: 99 994 813) and the STAG3 gene’s promoter region. Moreover, rs10241492 was identified as an eQTL for STAG3 from CommonMind in brain tissues [47]. In addition, STAG3 was predicted to be the target gene for rs10241492 from a previous GWAS study of EDU [58]. These results together showed that the STAG3 gene is potentially a target gene for rs10241492, specifically in microglia, consistent with findings in the original paper by Lee et al. [58].

Last but not least, Fig. S11 illustrates enhancer–promoter interactions identified specifically in microglia and neurons on chromosome 10. One microglia-specific enhancer–promoter interaction links a SCZ-associated GWAS locus in an active enhancer to SFXN2’s promoter. At the same time, SFXN2 is highly expressed only in microglia (FPKM = 4.41 in microglia versus ≤ 1 in the other three cell types). While a previous study suggested that SFXN2 was a potential SCZ risk gene because of linkage or pleiotropic effects [70], our results further indicated that SFXN2 is the putative target gene for this particular SCZ-associated GWAS locus. Additionally, multiple neuron-specific enhancer–promoter interactions connect SCZ-associated GWAS loci to INA (FPKM = 46.21 in neurons versus ≤ 1 in the other three cell types), suggesting that this gene is a potential novel target gene for this locus with the evidence that corresponding GWAS variants (rs11191557, rs11191558, rs11191559, rs10883832, and rs12413046) in this locus were also CommonMind eQTLs for INA.

Taken together, these examples showcase how SnapHiC-G identified enhancer–promoter interactions can aid the interpretation of noncoding GWAS variants and reveal underlying mechanistic insights. By integrating cell-type-specific gene expression data and epigenetic annotations, SnapHiC-G was able to decipher the critical roles of noncoding variants in disease etiology in relevant cell types.

Elucidating relevant cell types by heritability enrichment analysis

Next, we evaluated whether genetic heritability for complex diseases and traits was enriched for SNPs within the anchors of SnapHiC-G identified interactions in specific cell types. Using the same set of GWAS summary statistics data together with two other complex traits: body mass index (BMI) [71] and white blood cell count (WBC) [72], we performed stratified linkage disequilibrium score regression (S-LDSC) analysis [73]. In brief, S-LDSC estimates the proportion of SNP heritability from predefined SNP-level functional annotations using GWAS summary statistics, while accounting for linkage disequilibrium (LD) to identify functional categories enriched in SNP heritability and hence of functional relevance to the trait. In our study, the functional categories correspond to enhancer–promoter interactions from SnapHiC-G results in brain cell types, and our goal is to identify disease-relevant cell types for GWAS traits.

We first obtained SnapHiC-G identified interactions from astrocytes, microglia, oligodendrocytes, and L2/3 neurons, all with 261 cells, to construct functional categories. As previously described, SnapHiC-G requires at least one bin overlapping with a TSS to ensure that the algorithm captures promoter-anchored interactions. Therefore, both bins may overlap with a TSS. However, to distinguish bins that overlap with a TSS and those that do not overlap with a TSS, we focused on the case where only one bin overlaps with a TSS. We selected significant interactions in the four brain cell types to partition the genome. We defined SnapHiC-G anchor regions as the bins that overlap with TSS and SnapHiC-G target regions as the bins that do not overlap with TSS for each cell type. GWAS SNPs were annotated based on whether they fall into SnapHiC-G target regions for each cell type. We evaluated whether SNPs located in SnapHiC-G target regions show enriched heritability, in specific cell types. We similarly assessed heritability enrichment based on SnapHiC interactions. Figure 6 and Fig. S12 show SNP heritability enrichment results for GWAS variants from the 12 traits analyzed, using SnapHiC-G (Fig. 6) and SnapHiC (Fig. S12) target bins in four brain cell types. In most situations, SnapHiC-G demonstrates more significant P-values compared to SnapHiC, while SnapHiC exhibits larger enrichment scores than SnapHiC-G. To adjust for multiple testing, we used the Bonferroni-corrected P-value threshold of 0.05/48 = 0.001. Overall, we observed remarkably similar patterns of SNP heritability enrichment using SnapHiC-G or SnapHiC. First, AD SNP heritability was most strongly and significantly enriched in microglia-marked regions, which is consistent with the rich literature supporting vital roles played by microglia in the pathogenesis of AD [74]. The trends of magnitude for enrichment scores using SnapHiC-G and SnapHiC results were the same for AD. However, SnapHiC-G showed significant results (P = 1.07 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\times$\end{document} 10–6) for AD SNP heritability in microglia, while SnapHiC showed nonsignificant results (P = 0.01). Second, we observed that two neuropsychiatric traits with high genetic correlation, namely, BP and SCZ, have enrichment scores and significance in the same order, where regions marked in neurons were most strongly enriched, followed by astrocytes and oligodendrocytes. In contrast, enrichment in microglia was much lower. Third, our findings revealed strong enrichment in oligodendrocytes for PD, consistent with the literature [75]. Finally, results for white blood cell count showed that microglia were most strongly enriched, in agreement with the crucial roles that white blood cells have in the immune system.

Figure 6 LDSC heritability enrichment analysis results. Cells were down-sampled to match the number of cells in L2/3 neurons. The interactions were identified from 261 cells from each cell type. Numbers in the figure represent the enrichment score and colors represent the significance level in the −log10(P-value) scale. A higher enrichment score represents stronger enrichment in the corresponding cell type. Abbreviations: Alzheimer’s disease (AD), attention-deficit/hyperactivity disorder (ADHD), autism spectrum disorder (ASD), schizophrenia (SCZ), bipolar disorder (BP), educational attainment (EDU), intelligent quotient (IQ), body mass index (BMI), Parkinson’s disease (PD), major depressive disorder (MDD), neuroticism (NEU), white blood cell count (WBC). Interactions were identified using SnapHiC-G.

Discussion and Conclusion

In this work, we developed SnapHiC-G, a computational pipeline to detect enhancer–promoter interactions from scHi-C data based on the global background model. Our previous work, SnapHiC, the first computational pipeline to detect chromatin loops in scHi-C data, utilizes a combination of global and local background models to identify loop summits, therefore having limited power to detect enhancer–promoter interactions. More details about the difference between SnapHiC-G and SnapHiC are discussed in the Supplementary Text. Other existing global background methods, such as FitHiC2, HiC-DC+, FastHiC, Chromosight, and HiC-ACT, were designed for bulk Hi-C data and lacked the power to analyze sparse scHi-C contact data. To overcome the data sparsity challenge in scHi-C data, SnapHiC-G first applies the RWR algorithm to the observed raw contacts and then constructs a normalized contact probability matrix against the 1D genomic distance. For each candidate bin pair, SnapHiC-G applies a one-sample t-test to test whether its average normalized contact frequency across single cells is significantly greater than zero (global background). Combined with epigenetic annotations such as enhancer and promoter marks, SnapHiC-G enables the profiling of cell-type-specific enhancer–promoter interactions by analyzing scHi-C data from multiple cell types.

Data from mESCs and human brain cortical cells showed that SnapHiC-G identified more significant interactions than existing global background-based methods designed for bulk Hi-C data with notably higher sensitivity. This was mainly due to the RWR imputation step, which significantly reduced the sparsity of the scHi-C data, especially when the number of cells was low. Considering that Higashi, Fast-Higashi [23, 24], and scVI-3D [25] are not applicable to impute the 10Kb resolution scHi-C data, we have chosen not to compare with these imputation methods in our analysis. However, users may consider alternative imputation methods in the future, instead of the RWR imputation we utilized. Aggregating a small number of single cells can result in very sparse pseudo-bulk Hi-C data, which may lead to poor performance for methods designed for bulk Hi-C. In general, the extremely high sparsity in scHi-C data poses great challenges for loop callers designed for bulk Hi-C data, leading to increased false positives and reduced sensitivity in detecting enhancer–promoter interactions. The sensitivity of several bulk Hi-C methods, including HiC-DC+, HiC-ACT, and FitHiC2, suffers greatly when the number of cells is small (e.g. 100 mESCs) (Table 1). With the default parameter setting, Chromosight demonstrated a higher sensitivity with 100 mESCs than those methods and its performance was comparable to 742 mESCs. The primary distinction of Chromosight from other methods lies in its integration of concepts from computer vision to leverage visual patterns in Hi-C contact maps. This sets it apart from traditional statistical, probabilistic, or regression-based approaches.

We note that the number of interactions called by SnapHiC-G, although seemingly large and warrant further investigation, was comparable to existing methods for bulk data [28, 76, 77]. We further evaluated the precision of top interactions identified by each method against a reference list of interactions and showed that SnapHiC-G did not suffer from a higher FDR compared with other methods, demonstrating the usefulness of SnapHiC-G to reveal many missed enhancer–promoter interactions. Moreover, using experimentally tested enhancer–promoter connections in mESCs as true negatives, SnapHiC-G demonstrated a reasonable specificity of 81.5%. We acknowledge that the performance evaluation of SnapHiC-G and other methods depends on the reference loop list, where both true positives and true negatives are not gold standard and challenging to define. Moreover, the underlying chromatin interaction state might not be strictly dichotomous, posing further challenges in benchmarking loop callers on the genome-wide scale.

Although SnapHiC-G is designed to identify enhancer–promoter interactions, we found that it can identify a small proportion of structural loops using the 742 mESC scHi-C data (Fig. S4; Table S1). However, users are recommended to run other loop callers, such as SnapHiC, if their key interest is to identify structural loops instead of enhancer–promoter interactions.

Using GWAS summary statistics, we have also demonstrated the utility of enhancer–promoter interactions identified from SnapHiC-G in human brain cell types for prioritizing functionally relevant genes and cell types for complex human traits and diseases. SnapHiC-G can predict putative target genes for GWAS variants and help elucidate functionally relevant cell types in complex diseases.

For the human brain cell types, SnapHiC-G did not show substantial improvement over other methods compared with mESCs, especially when examining the most significant interactions for oligodendrocytes (Fig. S2; Fig. S3). Several reasons can explain this. First, there was no gold-standard data for evaluating cell-type-specific enhancer–promoter interactions. The reference interaction lists inferred from H3K4me3 PLAC-seq data were suggestive rather than optimal, meaning that the reference might miss many cell-type-specific enhancer–promoter interactions. Second, SnapHiC-G identified interactions had much significant FDRs than other bulk Hi-C methods. For identified oligodendrocyte interactions, the median of -log10(FDR) was 17.5 for SnapHiC-G, while the medians of other bulk Hi-C methods ranged from 1.7 to 7.5. The highly significant results of SnapHiC-G were due to both the RWR imputation and its global background nature; therefore, it can be hard to distinguish among identified interactions by ranking them by significance. Third, advantages of SnapHiC-G over methods developed for bulk Hi-C data were not particularly evident due to the relatively large (>1000) number of cells for oligodendrocytes. After aggregating the single cells to construct bulk Hi-C data, we obtained a comparable sequencing depth to traditional bulk Hi-C data with ~278 million intra-chromosomal contacts>20Kb. Consistent with SnapHiC, we observed a more significant power gain from SnapHiC-G when the number of cells was relatively small.

SnapHiC-G performs the statistical test across all cells from the same cell type, which means that signals from the input cells are aggregated, and the identified enhancer–promoter interactions are still at the population level, similar to bulk Hi-C analysis. However, chromatin folding can be highly variable and dynamic even among cells of similar identities [18], which is an exciting future direction. Enhancer–promoter contacts can be both cell-type-specific and developmental-stage-specific, suggesting that chromatin re-organization, coordinated with the cell cycle, plays a key role in shaping the 3D genome structure [78]. As more scHi-C data become available, multimodal data integration becomes another promising area of research to study the complexity and heterogeneity of chromatin interactions in single cells and can provide new insights into the regulatory mechanisms that underlie gene expression and cell differentiation [79]. For example, recent advances in single cell multimodal technologies [78, 80, 81] enable the simultaneous profiling of gene expression and chromatin spatial organization in the same single cell. The analysis of such multimodal single-cell datasets revealed the widespread chromatin rewiring before transcription activation, indicating that specific chromatin interactions are closely associated with transcriptional regulation and cell function during lineage specification. scHi-C data also have great potential for predicting structural variations in cancer genomes, which is beyond the scope of this work [82].

In our analysis, cell-type-specific epigenetic data considerably narrowed down candidate bin pairs. While such data aid the detection of cell-type-specific enhancer–promoter interactions, when they are not available, users can input only the TSS files to define the promoter regions and apply SnapHiC-G to identify promoter-interacting regions. With the rapid development of new technologies and more data available, scHi-C can trigger the study of fundamental questions about chromatin spatial organizations in individual cells during development, cancer cells, and different organs. Being able to identify cell-type-specific enhancer–promoter interactions from scHi-C data, SnapHiC-G results can be combined with the widely available GWAS results and epigenetic data, and it has a great potential to facilitate the discovery of regulatory chromatin interactions that are important for gene regulation in biologically relevant cell types.

Materials and Methods

SnapHiC-G algorithm

Step A. Imputation of contact probability using random walk with restart (RWR) with a sliding window approximation

We followed SnapHiC for imputing intrachromosomal contact probability in every cell using the RWR algorithm following scHiCluster [22]. Each autosome was divided into consecutive 10Kb bins, and each bin pair was converted to a binary representation, with 1 representing a nonzero contact and 0 representing no contact observed. An unweighted and undirected graph models each chromosome by defining bins as the nodes and adjacent bin pairs or bin pairs with a nonzero contact as edges. The RWR algorithm was then employed with a restart probability of 0.05 to estimate the likelihood of traveling between two nodes, allowing for the imputation of contact probability between all intrachromosomal bin pairs. The random-walk step captures information from global network structures, while the restart step captures information from local network structures. Following SnapHiC2, we adopted a sliding window approach when imputing missing contacts to reduce computational costs. Specifically, instead of performing RWR over the entire chromosome, we divided the original contact matrix into partially overlapping matrices of size 2Mb by 2Mb along the diagonal line with overlapping areas of 1Mb by 1Mb. We then performed RWR for all 10Kb bin pairs within each 2Mb by 2Mb submatrix along the diagonal to approximate contact probability. Only imputed contact probability in the middle rectangle areas was kept to avoid artifacts near corners.

Step B. Normalization of contact probability using 1D genomic distance

All bin pairs were stratified based on the genomic distance between two bins to account for the dependency between imputed contact probability and 1D genomic distance. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{k}=1,2,\dots, \mathrm{K}$\end{document} be the index of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{K}$\end{document} input cells. For a bin pair (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}) in cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{k}$\end{document} with a genomic distance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{d}$\end{document}, let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{A}}_{\mathrm{d}}^{\left(\mathrm{k}\right)}$\end{document} represent the strata including bin pairs in cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{k}$\end{document} with 1D genomic distance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{d}$\end{document}. The normalized contact probability (Z-score) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{z}}_{\mathrm{ij}}^{\left(\mathrm{k}\right)}$\end{document} was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{z}}_{\mathrm{ij}}^{\left(\mathrm{k}\right)}=\left({\mathrm{x}}_{\mathrm{ij}}^{\left(\mathrm{k}\right)}-{\mathrm{\mu}}_{\mathrm{d}}^{\left(\mathrm{k}\right)}\right)/{\mathrm{\sigma}}_{\mathrm{d}}^{\left(\mathrm{k}\right)}$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{x}}_{\mathrm{ij}}^{\left(\mathrm{k}\right)}$\end{document} is the contact probability between bin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i}$\end{document} and bin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{j}$\end{document} in cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{k}$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{\mu}}_{\mathrm{d}}^{\left(\mathrm{k}\right)}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{\sigma}}_{\mathrm{d}}^{\left(\mathrm{k}\right)}$\end{document} are mean and standard deviation of the contact probability of bin pairs in cell k within the strata \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{A}}_{\mathrm{d}}^{\left(\mathrm{k}\right)}$\end{document}.

Step C. Filtering for candidate chromatin interactions

First, we define the ‘AND’ bin pair as both sides overlapping with TSS and genes’ promoter regions, determined from the user’s input file or +/−500 bp of TSS. Next, we define the ‘XOR’ bin pair as only one side overlapping with TSS and genes’ promoter regions while the other side overlaps with enhancer regions. Chromatin interactions categorized as ‘AND’ or ‘XOR’ are the candidate SnapHiC-G enhancer–promoter interactions. In our analysis, we used TSS regions and H3K4me3 ChIP-seq or ATAC-seq peaks to define promoters and H3K27ac ChIP-seq peaks or ATAC-seq peaks to define putative enhancers [60, 83]. When these data are not available, users may provide only TSS annotations to define promoter regions and apply SnapHiC-G to identify promoter-interacting regions. Furthermore, candidate enhancer–promoter interactions with low mappability score (≤0.8) or overlapping with the ENCODE blacklist regions (mm10: http://mitra.stanford.edu/kundaje/akundaje/release/blacklists/mm10-mouse/mm10.blacklist.bed.gz; hg19: https://www.encodeproject.org/files/ENCFF001TDO/) were excluded. Our prior work was used to calculate each 10Kb bin’s sequence mappability [84].

Step D. Detecting enhancer–promoter interaction candidates

For each bin pair (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}), we applied the one-sample t-test across all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{K}$\end{document} cells to evaluate whether the contact probability was significantly higher than zero. We converted one sample t-test P-values into FDRs, again stratified by 1D genomic distance. We defined a bin pair as a significant chromatin interaction if (i) the average of z-scores across all input cells>0, (ii) the proportion of cells with a Z-score>1.96 (outlier cells)>10%, (iii) FDR<10%, and (iv) t-test statistic>3.

The computational cost of SnapHiC-G

Adopting a similar parallel computing strategy, SnapHiC-G is highly efficient regarding memory and computational time. We evaluated the computational cost of SnapHiC-G-specific steps (Step C and Step D) with both mouse and human scHi-C data under different settings and summarized the results in Table S5.

Processing scHi-C data

We used the same procedure as in SnapHiC to process the single-cell Hi-C data of mESCs and the human prefrontal cortex sn-m3C-seq data. For mESCs, we aligned the raw read pairs in fastq format to the mm10 genome, removed duplications, and then chose the top 742 cells (>150 000 contacts in each cell) for downstream analysis. For the sn-m3C-seq data from human cortex, we used reference genome hg19 to process the data and then removed duplications. We chose the top 2869 cells (>150 000 contacts in each cell) for downstream analysis. We only analyzed cells from the same cell type together and used cell type annotations reported by the original study [45].

Down-sampling of scHi-C data

We randomly permuted 742 quality-controlled cells for the mESC scHi-C data and then performed down-sampling by selecting the first 100 and 300 cells from the pool of 742 cells. As for the astrocytes (338 cells), microglia (323 cells), and oligodendrocytes (1038 cells) from the human prefrontal cortex sn-m3C-seq data, we permuted the cells and selected the first 261 cells from each cell type to match the number of L2/3 neurons.

Identification of SnapHiC-G chromatin enhancer–promoter interactions

SnapHiC-G was applied to 742 mESCs (and downsampled to 100 and 300 mESCs) scHi-C data and four cell types of the human prefrontal cortex sn-m3C-seq data to identify 10Kb bin chromatin enhancer–promoter interactions on autosomal chromosomes with 20Kb to 1Mb distance range. Bin pairs within 20Kb were excluded from analyses. Two major reasons support SnapHiC-G to search for long-range enhancer–promoter interactions beyond 20Kb. First, enhancer and promoter in short range (<20Kb) may interact with each other due to random collision, instead of chromatin looping. Therefore, Hi-C data are not informative to detect such short-range (<20Kb) interactions. More importantly, many studies have demonstrated that cis-regulatory elements control the expression of distal genes via long-range chromatin interactions [2, 85], and Hi-C has been widely used to measure such long-range looping between enhancer and promoter. Therefore, we only report enhancer–promoter interactions in the range of 20Kb to 1Mb in this study.

Existing methods for bulk Hi-C data

Many methods, including FitHiC2, FastHiC, HiC-ACT, HiC-DC+, and Chromosight, have been developed for identifying long-range chromatin interactions from bulk Hi-C data. Specifically, FitHiC2 [28, 86] is a spline regression-based approach to identify intra-chromosomal chromatin interactions. FastHiC [29] is a Bayesian hidden Markov random field method that models the spatial dependency structure in high-resolution Hi-C data for detecting biologically meaningful chromatin interactions (i.e. peak calling). HiC-ACT [30] is an aggregated Cauchy test–based approach that combines P-values from other methods (e.g. FitHiC2) without knowing the underlying spatial dependency structure. HiC-DC+ [31] is a negative binomial regression–based method based on genomic distance, guanine-cytosine (GC) content, and mappability features to predict chromatin interactions. Chromosight treats a Hi-C contact map as an image (with pixels corresponding to pairwise contacts) and employs a convolution algorithm from the computer vision domain to detect loops [27].

Identification of loops/interactions using other Hi-C methods

We identified chromatin loops/interactions with HiCCUPS, Chromosight, FastHiC, FitHiC2, HiC-ACT, and HiC-DC+ from aggregated contact matrix and SnapHiC from single-cell contact matrices. To apply the methods developed for bulk Hi-C data, we generated pseudo-bulk Hi-C data by aggregating 10Kb-resolution scHi-C contact matrices across single cells (i.e., sum up the contacts) [32]. Then, we applied Chromosight, FitHiC2, FastHiC, HiC-ACT, and HiC-DC+ to the pseudo-bulk Hi-C data with lenient significance thresholds to account for the sparsity of scHi-C data. We used the following criteria for bulk Hi-C methods to identify significant interactions: FDR<10% from FitHiC2; posterior probability>0.9 from FastHiC; local neighborhood smoothed P-values<10−6 from HiC-ACT; FDR< 10% from HiC-DC+; and Pearson correlation coefficient >0.3 from Chromosight. For SnapHiC, we used FDR<10% and t-statistics>3 to select significant interactions. We further required interactions to be within the 20Kb–1Mb 1D genomic distance, with high mappability (>0.8) and no overlapping with ENCODE blacklist regions for standard quality control. To fairly compare with SnapHiC-G-identified interactions, we required interactions to be ‘AND’ or ‘XOR’ as an additional filtering criterion. Note that the number of interactions for SnapHiC was reduced compared with the original SnapHiC paper because of this additional filtering criterion.

Definition of overlapped enhancer–promoter interactions

We followed the same definition of overlapped loops as in SnapHiC. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{d}}_{\mathrm{im}}$\end{document} denote the 1D genomic distance between the center of bin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i}$\end{document} and the center of bin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{m}$\end{document}, and we define the distance between (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}) and (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{m},\mathrm{n}$\end{document}) as the maximum of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{d}}_{\mathrm{im}}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{d}}_{\mathrm{jn}}$\end{document}. For an enhancer–promoter interaction (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}), if there exists an enhancer–promoter interaction (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{m},\mathrm{n}$\end{document}) in set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{S}$\end{document} such that the distance between (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}) and (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{m},\mathrm{n}$\end{document}) is within 20Kb, we define that the enhancer–promoter interaction (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}) overlaps with the set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{S}$\end{document}.

Definition of exact overlapped enhancer–promoter interactions

For an enhancer–promoter interaction (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}), if there exists an enhancer–promoter interaction (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{m},\mathrm{n}$\end{document}) in set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{S}$\end{document} such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i}=\mathrm{m}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{j}=\mathrm{n}$\end{document}, we define that the enhancer–promoter interaction (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{i},\mathrm{j}$\end{document}) exactly overlaps with the set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{S}$\end{document}.

Definition of reference interaction lists for evaluation

The reference interactions for mESCs consist of 10Kb interactions called from HiCCUPS with deeply sequenced bulk Hi-C data [35] and MAPS-identified significant interactions from H3K4me3 PLAC-seq [36], cohesin HiChIP [37], and H3K27ac HiChIP data [38]. The reference interactions for oligodendrocytes, microglia, and neurons from the human prefrontal cortex were constructed using MAPS-identified H3K4me3 PLAC-seq data from corresponding cell types [60], following SnapHiC. The same filtering steps applied to SnapHiC-G and other methods were applied to the reference lists to ensure a fair evaluation.

Definition of cell-type-specific SnapHiC-G enhancer–promoter interactions

Cell-type-specific enhancer–promoter interactions were a subset of the SnapHiC-G enhancer–promoter interactions detected from (downsampled) 261 cells in each of the four cell types: astrocytes, oligodendrocytes, microglia, and L2/3 excitatory neurons. Specifically, an enhancer–promoter interaction detected from a cell type was defined as a cell-type-specific SnapHiC-G enhancer–promoter interaction if none of the other three cell types’ identified enhancer–promoter interactions overlapped with it (Definition of overlapped enhancer–promoter interactions).

Enrichment of eQTL–TSS pairs

We evaluated whether SnapHiC-G-identified interactions from the human cortex data are enriched in eQTL-TSS bin pairs, specifically “brain cell type” eQTL-TSS bin pairs [46], CommonMind brain eQTL-TSS bin pairs [47], and GTEx consortium liver eQTL-TSS bin pairs [48], for all four cell types that we considered. Specifically, for eQTL-TSS bin pairs in brain cell types, we used eQTLs with Bonferroni-corrected P-values<0.05 within each cell type. Because the brain cell type eQTL data did not have results in L2/3 neurons, we used eQTL results of excitatory neurons to compare with SnapHiC-G identified bin pairs in L2/3 neurons. For each SnapHiC-G identified interaction, we constructed a matched pseudo bin pair as a control: if only one bin contains the promoter, we kept the bin with the promoter as the center and flipped the other bin to be on the opposite side of the center but with the same distance from the center; if both bins contain promoters, we randomly selected one bin with probability 0.5 and kept it as the center, and then flipped the other bin to the other side of the center. Next, we removed the duplicates between controls and SnapHiC-G identified bin pairs.

We constructed a two-by-two table for the union of SnapHiC-G identified interactions and pseudo bin pairs, categorizing each bin pair by whether it is a SnapHiC-G identified interaction or an eQTL-TSS bin pair. The P-value for independence between these two features was calculated using a two-sided Fisher’s exact test.

Gene expression analysis at cell-type-specific enhancer–promoter interactions

The FPKM values of each protein-coding gene in human astrocytes, neurons, microglia, and oligodendrocytes were acquired from Zhang et al. [49]. We used average FPKM values across biological replicates of the same cell type to quantify cell-type-specific gene expression levels.

Selection of GWAS SNPs associated with neuropsychiatric disorders and complex traits

First, we gathered significant (P-values< 5 × 10−8) GWAS SNPs for 10 brain-related traits, including AD [50], ADHD [51], ASDs [52], bipolar disorder (BP) [53], SCZ [54], PD [55], MDD [56], NEU [57], EDU [58], and IQ [59]. Next, we overlapped these GWAS SNPs with active enhancers in four cell types (astrocytes, neurons, microglia, and oligodendrocytes), resulting in 9764 SNP–trait associations (8516 unique GWAS SNPs).

Web resources for figure generation

Figures 1, 2, 4, 6, S2–S8, and S12 were created with the R package ggplot2 (https://CRAN.R-project.org/package=ggplot2). Figure 3A was created with the R package forestplot (https://CRAN.R-project.org/package=forestplot). Figure 3B was created with the R package UpSetR (https://CRAN.R-project.org/package=UpSetR). Figures 5, S10, and S11 were created with IGV (https://igv.org/app). Figures 1, 5, S2, S10, and S11 were arranged by BioRender (https://biorender.com/). Figure S9 was created with Juicebox (https://aidenlab.org/juicebox/).

Key Points

We propose a new method, SnapHiC-G, to identify long-range enhancer–promoter interactions from single-cell Hi-C data.

SnapHiC-G overcomes the limited sensitivity of previously developed methods.

SnapHiC-G can significantly enhance the interpretation and prioritization of GWAS variants.

Supplementary Material

bib_supp_08092024_bbae426

Acknowledgements

We thank assistance from the UNC Intellectual and Developmental Disabilities Research Center (NICHD; P50 HD103573).

 

Conflict of interest: W.Z. is an employee at Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc., Rahway, NJ, USA.

Funding

This work was supported by the National Institutes of Health [R35HG011922 to M.H., U01DA052713 to Y.L., R01MH125236 to Y.L., P50HD103573 to Y.L., and R01NR019245 to W.L.].

Data availability

The proposed SnapHiC-G method is implemented in the SnapHiC-G Python package, freely available on GitHub: https://github.com/wjzhong/SnapHiC-G.
==== Refs
References

1. Fulco  CP, Munschauer  M, Anyoha  R. et al.  Systematic mapping of functional enhancer–promoter connections with CRISPR interference. Science  2016;354 :769–73.27708057
2. Zhong  W, Liu  W, Chen  J. et al.  Understanding the function of regulatory DNA interactions in the interpretation of non-coding GWAS variants. Front Cell Dev Biol  2022;10 :957292.36060805
3. Li  Y, Hu  M, Shen  Y. Gene regulation in the 3D genome. Hum Mol Genet  2018;27 :R228–33.29767704
4. Fudenberg  G, Imakaev  M, Lu  C. et al.  Formation of chromosomal domains by loop extrusion. Cell Rep  2016;15 :2038–49.27210764
5. Rowland  B, Huh  R, Hou  Z. et al.  THUNDER: a reference-free deconvolution method to infer cell type proportions from bulk hi-C data. PLoS Genet  2022;18 :e1010102.35259165
6. Foster  HA, Bridger  JM. The genome and the nucleus: a marriage made by evolution. Genome organisation and nuclear architecture. Chromosoma  2005;114 :212–29.16133352
7. Nagano  T, Lubling  Y, Stevens  TJ. et al.  Single-cell hi-C reveals cell-to-cell variability in chromosome structure. Nature  2013;502 :59–64.24067610
8. Cattoni  DI, Cardozo Gizzi  AM, Georgieva  M. et al.  Single-cell absolute contact probability detection reveals chromosomes are organized by multiple low-frequency yet specific interactions. Nat Commun  2017;8 :1753.29170434
9. Finn  EH, Pegoraro  G, Brandão  HB. et al.  Extensive heterogeneity and intrinsic variation in spatial genome organization. Cell  2019;176 :1502–1515.e10.30799036
10. Bintu  B, Mateo  LJ, Su  J-H. et al.  Super-resolution chromatin tracing reveals domains and cooperative interactions in single cells. Science  2018;362 :eaau1783.30361340
11. Collombet  S, Ranisavljevic  N, Nagano  T. et al.  Parental-to-embryo switch of chromosome organization in early embryogenesis. Nature  2020;580 :142–6.32238933
12. Wang  S, Su  JH, Beliveau  BJ. et al.  Spatial organization of chromatin domains and compartments in single chromosomes. Science  2016;353 :598–602.27445307
13. Boettiger  AN, Bintu  B, Moffitt  JR. et al.  Super-resolution imaging reveals distinct chromatin folding for different epigenetic states. Nature  2016;529 :418–22.26760202
14. Finn  EH, Misteli  T. Molecular basis and biological function of variability in spatial genome organization. Science  2019;365 :eaaw9498.31488662
15. Galitsyna  AA, Gelfand  MS. Single-cell hi-C data analysis: safety in numbers. Brief Bioinform  2021;22 :bbab316.34406348
16. Zhou  T, Zhang  R, Ma  J. The 3D genome structure of single cells. Annu Rev Biomed data Sci  2021;4 :21–41.34465168
17. Ramani  V, Deng  X, Qiu  R. et al.  Massively multiplex single-cell hi-C. Nat Methods  2017;14 :263–6.28135255
18. Nagano  T, Lubling  Y, Várnai  C. et al.  Cell-cycle dynamics of chromosomal organization at single-cell resolution. Nature  2017;547 :61–7.28682332
19. Tan  L, Xing  D, Chang  C-H. et al.  Three-dimensional genome structures of single diploid human cells. Science  2018;361 :924–8.30166492
20. Nguyen  HQ, Chattoraj  S, Castillo  D. et al.  3D mapping and accelerated super-resolution imaging of the human genome using in situ sequencing. Nat Methods  2020;17 :822–32.32719531
21. Ramani  V, Deng  X, Qiu  R. et al.  Sci-hi-C: a single-cell hi-C method for mapping 3D genome organization in large number of single cells. Methods  2020;170 :61–8.31536770
22. Zhou  J, Ma  J, Chen  Y. et al.  Robust single-cell hi-C clustering by convolution- and random-walk-based imputation. Proc Natl Acad Sci U S A  2019;116 :14011–8.31235599
23. Zhang  R, Zhou  T, Ma  J. Multiscale and integrative single-cell hi-C analysis with higashi. Nat Biotechnol  2022;40 :254–61.34635838
24. Zhang  R, Zhou  T, Ma  J. Ultrafast and interpretable single-cell 3D genome analysis with fast-higashi. Cell Syst  2022;13 :798–807.36265466
25. Zheng  Y, Shen  S, Keleş  S. Normalization and de-noising of single-cell hi-C data with BandNorm and scVI-3D. Genome Biol  2022;23 :222.36253828
26. Liu  T, Wang  Z. scHiCEmbed: bin-specific embeddings of single-cell hi-C data using graph auto-encoders. Genes (Basel)  2022;13 :1048.35741810
27. Matthey-Doret  C, Baudry  L, Breuer  A. et al.  Computer vision for pattern detection in chromosome contact maps. Nat Commun  2020;11 :5795.33199682
28. Kaul  A, Bhattacharyya  S, Ay  F. Identifying statistically significant chromatin contacts from hi-C data with FitHiC2. Nat Protoc  2020;15 :991–1012.31980751
29. Xu  Z, Zhang  G, Wu  C. et al.  FastHiC: a fast and accurate algorithm to detect long-range chromosomal interactions from hi-C data. Bioinformatics  2016;32 :2692–5.27153668
30. Lagler  TM, Abnousi  A, Hu  M. et al.  HiC-ACT: improved detection of chromatin interactions from hi-C data via aggregated Cauchy test. Am J Hum Genet  2021;108 :257–68.33545029
31. Sahin  M, Wong  W, Zhan  Y. et al.  HiC-DC+ enables systematic 3D interaction calls and differential analysis for hi-C and HiChIP. Nat Commun  2021;12 :3366.34099725
32. Yu  M, Abnousi  A, Zhang  Y. et al.  SnapHiC: a computational pipeline to identify chromatin loops from single-cell hi-C data. Nat Methods  2021;18 :1056–9.34446921
33. Li  X, Lee  L, Abnousi  A. et al.  SnapHiC2: a computationally efficient loop caller for single cell hi-C data. Comput Struct Biotechnol J  2022;20 :2778–83.35685374
34. Liu  W, Yang  Y, Abnousi  A. et al.  MUNIn: a statistical framework for identifying long-range chromatin interactions from multiple samples. HGG Adv  2021;2 :100036.34485947
35. Bonev  B, Mendelson Cohen  N, Szabo  Q. et al.  Multiscale 3D genome rewiring during mouse neural development. Cell  2017;171 :557–572.e24.29053968
36. Juric  I, Yu  M, Abnousi  A. et al.  MAPS: model-based analysis of long-range chromatin interactions from PLAC-seq and HiChIP experiments. PLoS Comput Biol  2019;15 :e1006982.30986246
37. Mumbach  MR, Rubin  AJ, Flynn  RA. et al.  HiChIP: efficient and sensitive analysis of protein-directed genome architecture. Nat Methods  2016;13 :919–22.27643841
38. Mumbach  MR, Satpathy  AT, Boyle  EA. et al.  Enhancer connectome in primary human cells identifies target genes of disease-associated DNA elements. Nat Genet  2017;49 :1602–12.28945252
39. Zhou  HY, Katsman  Y, Dhaliwal  NK. et al.  A Sox2 distal enhancer cluster regulates embryonic stem cell differentiation potential. Genes Dev  2014;28 :2699–711.25512558
40. Li  Y, Rivera  CM, Ishii  H. et al.  CRISPR reveals a distal super-enhancer required for Sox2 expression in mouse embryonic stem cells. PloS One  2014;9 :e114485.25486255
41. Engreitz  JM, Haines  JE, Perez  EM. et al.  Local regulation of gene expression by lncRNA promoters, transcription and splicing. Nature  2016;539 :452–5.27783602
42. Moorthy  SD, Davidson  S, Shchuka  VM. et al.  Enhancers and super-enhancers have an equivalent regulatory role in embryonic stem cells through regulation of single or multiple genes. Genome Res  2017;27 :246–58.27895109
43. Fulco  CP, Nasser  J, Jones  TR. et al.  Activity-by-contact model of enhancer–promoter regulation from thousands of CRISPR perturbations. Nat Genet  2019;51 :1664–9.31784727
44. Ren  B . ENCSR000CCB. ENCODE Datasets  2011.
45. Lee  D-S, Luo  C, Zhou  J. et al.  Simultaneous profiling of 3D genome structure and DNA methylation in single human cells. Nat Methods  2019;16 :999–1006.31501549
46. Bryois  J, Calini  D, Macnair  W. et al.  Cell-type-specific cis-eQTLs in eight human brain cell types identify novel risk genes for psychiatric and neurological disorders. Nat Neurosci  2022;25 :1104–12.35915177
47. Fromer  M, Roussos  P, Sieberts  SK. et al.  Gene expression elucidates functional impact of polygenic risk for schizophrenia. Nat Neurosci  2016;19 :1442–53.27668389
48. GTEx Consortium . Genetic effects on gene expression across human tissues. Nature  2017;550 :204–13.29022597
49. Zhang  Y, Sloan  SA, Clarke  LE. et al.  Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse. Neuron  2016;89 :37–53.26687838
50. Schwartzentruber  J, Cooper  S, Liu  JZ. et al.  Genome-wide meta-analysis, fine-mapping and integrative prioritization implicate new Alzheimer’s disease risk genes. Nat Genet  2021;53 :392–402.33589840
51. Demontis  D, Walters  RK, Martin  J. et al.  Discovery of the first genome-wide significant risk loci for attention deficit/hyperactivity disorder. Nat Genet  2019;51 :63–75.30478444
52. Grove  J, Ripke  S, Als  TD. et al.  Identification of common genetic risk variants for autism spectrum disorder. Nat Genet  2019;51 :431–44.30804558
53. Stahl  EA, Breen  G, Forstner  AJ. et al.  Genome-wide association study identifies 30 loci associated with bipolar disorder. Nat Genet  2019;51 :793–803.31043756
54. Pardiñas  AF, Holmans  P, Pocklington  AJ. et al.  Common schizophrenia alleles are enriched in mutation-intolerant genes and in regions under strong background selection. Nat Genet  2018;50 :381–9.29483656
55. Pankratz  N, Beecham  GW, DeStefano  AL. et al.  Meta-analysis of Parkinson’s disease: identification of a novel locus, RIT2. Ann Neurol  2012;71 :370–84.22451204
56. Howard  DM, Adams  MJ, Clarke  T-K. et al.  Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions. Nat Neurosci  2019;22 :343–52.30718901
57. Nagel  M, Jansen  PR, Stringer  S. et al.  Meta-analysis of genome-wide association studies for neuroticism in 449,484 individuals identifies novel genetic loci and pathways. Nat Genet  2018;50 :920–7.29942085
58. Lee  JJ, Wedow  R, Okbay  A. et al.  Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals. Nat Genet  2018;50 :1112–21.30038396
59. Savage  JE, Jansen  PR, Stringer  S. et al.  Genome-wide association meta-analysis in 269,867 individuals identifies new genetic and functional links to intelligence. Nat Genet  2018;50 :912–9.29942086
60. Nott  A, Holtman  IR, Coufal  NG. et al.  Brain cell type-specific enhancer–promoter interactome maps and disease-risk association. Science  2019;366 :1134–9.31727856
61. Yang  X, Wen  J, Yang  H. et al.  Functional characterization of Alzheimer’s disease genetic variants in microglia. Nat Genet  2023;55 :1735–44.37735198
62. Trubetskoy  V, Pardiñas  AF, Qi  T. et al.  Mapping genomic loci implicates genes and synaptic biology in schizophrenia. Nature  2022;604 :502–8.35396580
63. Mollereau  C, Simons  MJ, Soularue  P. et al.  Structure, tissue distribution, and chromosomal localization of the prepronociceptin gene. Proc Natl Acad Sci U S A  1996;93 :8666–70.8710928
64. Darland  T, Heinricher  MM, Grandy  DK. Orphanin FQ/nociceptin: a role in pain and analgesia, but so much more. Trends Neurosci  1998;21 :215–21.9610886
65. Girgenti  MJ, Wang  J, Ji  D. et al.  Transcriptomic organization of the human brain in post-traumatic stress disorder. Nat Neurosci  2021;24 :24–33.33349712
66. Løkhammer  S, Stavrum  A-K, Polushina  T. et al.  An epigenetic association analysis of childhood trauma in psychosis reveals possible overlap with methylation changes associated with PTSD. Transl Psychiatry  2022;12 :177.35501310
67. Jordanovski  D, Herwartz  C, Pawlowski  A. et al.  The hypoxia-inducible transcription factor ZNF395 is controlled by IĸB kinase-signaling and activates genes involved in the innate immune response and cancer. PloS One  2013;8 :e74911.24086395
68. Sahu  A, Chowdhury  HA, Gaikwad  M. et al.  Integrative network analysis identifies differential regulation of neuroimmune system in schizophrenia and bipolar disorder. Brain, Behav Immun - Heal  2020;2 :100023.
69. Chen  W-T, Lu  A, Craessaerts  K. et al.  Spatial transcriptomics and In situ sequencing to study Alzheimer’s disease. Cell  2020;182 :976–991.e19.32702314
70. Wang  J-Y, Li  X-Y, Li  H-J. et al.  Integrative analyses followed by functional characterization reveal TMEM180 as a schizophrenia risk gene. Schizophr Bull  2021;47 :1364–74.33768244
71. Yengo  L, Sidorenko  J, Kemper  KE. et al.  Meta-analysis of genome-wide association studies for height and body mass index in ∼700000 individuals of European ancestry. Hum Mol Genet  2018;27 :3641–9.30124842
72. Vuckovic  D, Bao  EL, Akbari  P. et al.  The polygenic and monogenic basis of blood traits and diseases. Cell  2020;182 :1214–1231.e11.32888494
73. Finucane  HK, Bulik-Sullivan  B, Gusev  A. et al.  Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat Genet  2015;47 :1228–35.26414678
74. Hemonnot  A-L, Hua  J, Ulmann  L. et al.  Microglia in Alzheimer disease: well-known targets and new opportunities. Front Aging Neurosci  2019;11 :233.31543810
75. Bryois  J, Skene  NG, Hansen  TF. et al.  Genetic identification of cell types underlying brain complex traits yields insights into the etiology of Parkinson’s disease. Nat Genet  2020;52 :482–93.32341526
76. Bhattacharyya  S, Chandra  V, Vijayanand  P. et al.  Identification of significant chromatin contacts from HiChIP data by FitHiChIP. Nat Commun  2019;10 :4221.31530818
77. Li  X, An  Z, Zhang  Z. Comparison of computational methods for 3D genome analysis at single-cell hi-C level. Methods  2020;181–182 :52–61.
78. Liu  Z, Chen  Y, Xia  Q. et al.  Linking genome structures to functions by simultaneous single-cell hi-C and RNA-seq. Science  2023;380 :1070–6.37289875
79. Boninsegna  L, Yildirim  A, Polles  G. et al.  Integrative genome modeling platform reveals essentiality of rare contact events in 3D genome organizations. Nat Methods  2022;19 :938–49.35817938
80. Zhou  T, Zhang  R, Jia  D. et al.  Concurrent profiling of multiscale 3D genome organization and gene expression in single mammalian cells. bioRxiv Prepr Serv Biol  2023.
81. Wen  X, Luo  Z, Zhao  W. et al.  Single-cell multiplex chromatin and RNA interactions in ageing human brain. Nature  2024;628 :648–56.38538789
82. Wang  X, Luan  Y, Yue  F. EagleC: a deep-learning framework for detecting a full range of structural variations from bulk and single-cell contact maps. Sci Adv  2022;8 :eabn9215.35704579
83. Kubo  N, Ishii  H, Xiong  X. et al.  Promoter-proximal CTCF binding promotes distal enhancer-dependent gene activation. Nat Struct Mol Biol  2021;28 :152–61.33398174
84. Hu  M, Deng  K, Selvaraj  S. et al.  HiCNorm: removing biases in hi-C data via Poisson regression. Bioinformatics  2012;28 :3131–3.23023982
85. Liu  W, Zhong  W, Chen  J. et al.  Understanding regulatory mechanisms of brain function and disease through 3D genome organization. Genes (Basel)  2022;13 :586.35456393
86. Ay  F, Bailey  TL, Noble  WS. Statistical confidence estimation for hi-C data reveals regulatory chromatin contacts. Genome Res  2014;24 :999–1011.24501021
