Genetic variants shape single-cell chromatin accessibility and disease risk

Wait 5 sec.

MainUnderstanding how genetic variation influences complex traits and disease risk remains a major challenge in human genomics. Although genome-wide association studies (GWAS) have identified tens of thousands of risk loci, most lie in non-coding regions6, complicating interpretation. One common approach to bridge this gap is to map molecular QTLs (molQTLs), identifying genetic variants that affect intermediate molecular phenotypes such as gene expression or chromatin accessibility, and then test for colocalization with GWAS loci to infer potential regulatory mechanisms7. Most molQTL studies have focused on eQTLs derived from bulk RNA sequencing (RNA-seq), which have low sensitivity for cell-type-specific effects8,9. Recently, well-powered single-cell RNA-seq (scRNA-seq) studies have shown that genetic regulation of gene expression varies across cell types and states10,11,12,13,14. However, many GWAS signals still lack corresponding regulatory effects, a gap often termed the ‘missing regulation’ problem15,16,17.caQTLs offer a complementary modality that directly reflects the physical state of regulatory elements and captures non-coding variant effects on enhancer and promoter activity. Recent studies have suggested that caQTLs colocalize with more GWAS loci than eQTLs15,16,18,19, particularly at distal elements, but current caQTL maps remain limited in scale or resolution20,21, or rely on genotypes imputed from the chromatin data22. Many caQTLs also lack corresponding eQTLs for reasons that remain unclear16,23, and dynamic, cell-state-dependent regulatory effects are often overlooked11,24.Here we present a large-scale atlas of genetic effects on single-cell chromatin accessibility in peripheral blood mononuclear cells (PBMCs). We profiled over 4.6 million single-cell assay for transposase-accessible chromatin with sequencing (scATAC-seq) and 95,728 multiome (RNA + ATAC) nuclei across 1,042 individuals. After quality control, we analysed 3.5 million high-quality nuclei from 922 donors, with matched whole-genome sequencing data. We identified 243,225 caQTLs across 28 immune cell types, with more than half showing cell-type-specific effects. Integrating caQTL and eQTL data revealed over 70,000 colocalized signals, including 28,834 candidate cis-regulatory elements (cCREs) supported by Mendelian randomization. We further detected thousands of cell-state-dependent caQTLs and GWAS colocalizations missed by eQTLs alone, and found that combining caQTL and eQTL data improved gene regulatory networks (GRNs) inference over conventional multiome-based methods.Together with accompanying articles exploring the broader effects of common and rare genetic variation5 and tandem repeat variants25 on gene expression, this work provides a single-cell view of genetic regulation of chromatin accessibility in human immune cells, filling gaps left by eQTL-only approaches to reveal how non-coding variation shapes immune function and disease risk.The TenK10K multiome cohortWe generated the TenK10K phase I multiome dataset by profiling chromatin accessibility (scATAC-seq) in 4.6 million nuclei from 952 donors, as well as multiome (RNA + ATAC) from 90 donors (Methods and Supplementary Tables 1–6). We obtained the 952 donors from the TenK10K pilot cohort, known as OneK1K or the Tasmanian Ophthalmic Biobank10 (Methods). These data were complemented by 5.4 million PBMC scRNA-seq profiles from 1,925 donors5, including 932 that overlapped with matched scATAC-seq profiles (Methods and Extended Data Fig. 1a). Cell types were annotated via label transfer against a CITE-seq reference (Methods), yielding 28 immune cell types consistent with the scRNA-seq dataset (Extended Data Figs. 1b and 2 and Supplementary Tables 4–6). After quality controls, 3,472,552 scATAC-seq nuclei from 922 donors were retained for discovery analyses, with an additional 60 donors used for replication (Methods, Extended Data Fig. 1c and Supplementary Fig. 1). On average, each donor contributed 3,766 nuclei; CD14+ monocytes (CD14Mono) were the most abundant cell type, whereas proliferating CD8+ T cells (CD8Proliferating) were the rarest (Extended Data Fig. 1d).Chromatin peaks were called independently in each cell type and merged into a unified list of 440,996 peaks (Methods and Extended Data Fig. 1e). Peak widths ranged from 201 to 5,545 bp (median of 499 bp). Most peaks were distal to the promoter, with only 4.0% within 200 bp of a transcription start site (TSS; Extended Data Fig. 1f). Compared with existing chromatin accessibility datasets, only 44,015 peaks (9.98%) lacked overlap with known peaks or cCREs (Methods and Supplementary Table 7).To assess the biological specificity of these peaks, we clustered cell-type-enriched features and observed grouping by haematopoietic lineage (Extended Data Fig. 2c,d). Transcription factor motif enrichment analyses supported known lineage functions: for example, monocyte-specific peaks were enriched for CEBP family motifs (CEBPB and CEBPD), consistent with their established roles in monocyte differentiation and activation (Extended Data Fig. 2e). Together, these results confirm that chromatin accessibility profiles capture cell-type-specific regulatory signals across the immune landscape.Mapping common variant single-cell caQTLsWe mapped caQTLs across 28 immune cell types in 922 donors, testing 440,996 peaks against common variants (minor allele frequency (MAF) ≥ 5%; approximately 9 million SNPs and 3.4 million indels) within 1 Mb of the peak centre (Methods). Using permutation-based correction (q  0.8) involving 28,834 peaks and 11,204 genes (Fig. 2a and Supplementary Table 10). We define these peaks as cCREs, enabling peak–gene linkage without paired multiome data.Fig. 2: Characterization of candidate cCREs across immune cell types.Full size imagea, Cell-type-specific and shared cCREs identified by colocalizing caQTLs and eQTLs. The number of peak–gene pairs colocalized in or shared by multiple cell types (right). The number of colocalized peak–gene pairs (PPH4 > 0.8) per cell type, grouped by lineage in the same colour scheme (left). b, Combined violin and boxplots of the top variants for the caQTL (chr. 10: 45592479–45592785) and eQTL (MARCH8) in NK cells and cDC2, showing standardized pseudobulk chromatin accessibility and gene expression across the three genotype groups (n = 922 individuals). Slope estimates and P values from a two-sided Student’s t-test are annotated. c, Locus plot of the caQTL. The −log10(nominal P values) from a two-sided Student’s t-test against genomic position, with the two top variants circled in black and the peak as a blue vertical line (top). Gene track, with forward-strand genes in blue and reverse-strand genes in green (middle). Zoomed-in MARCH8 locus and chromatin peak map across five cell types (bottom). d, Zoomed-in gene and peak tracks. Each dash in the peak track is an open chromatin peak. The sky blue vertical line marks the caPeak (chr. 10: 45592479–45592785).Colocalization events varied by cell type, with classical monocytes (CD14Mono) showing the most events (12,940), including 8,402 CD14Mono-specific pairs and 1,743 shared with non-classical monocytes (CD16Mono) (Fig. 2a). For example, the inflammatory bowel disease (IBD) risk gene IRGM showed a regulatory effect via promoter accessibility (chr. 5: 150846177–150847636) in CD8TEM cells but not CD8+ naive T cells (CD8Naive; Supplementary Fig. 16). This cell-state-dependent pattern is consistent with the known role of IRGM as an IBD risk locus37 and illustrates how integrating caQTL and eQTL data can resolve context-specific regulatory activity.Although coloc infers whether two phenotypes share a causal variant, it does not estimate effect direction or magnitude. To address this, we applied summary-based Mendelian randomization (SMR)38, using chromatin accessibility as the exposure and gene expression as the outcome (Methods). Testing 3,667,924 peak–gene–cell-type combinations (matching those tested in coloc), we identified 81,550 significant associations at PSMR  0.8), caPeaks by colocalization with eQTL (PPH4 > 0.8) and set focal regions to ±1 Mb of the gene body. We applied this to IBD, one of the best-powered autoimmune GWAS, across all immune cell types.Across 21 immune cell types, we identified 1,085 credible set–phenotype pairs, of which 882 were single-variant credible sets, among the 916 matching focal gene–cell-type pairs (Supplementary Table 15 and Supplementary Fig. 23). For example, 11:1852991:C:T, 2 kb upstream of the IBD-associated gene LSP1, was fine-mapped as a causal variant for IBD risk. In CD14 monocytes, this variant was strongly associated with both LSP1 expression (PeQTL = 2.21 × 10−131) and accessibility at chr. 11: 1852420–1853362 (PcaQTL = 9.77 × 10−47; Fig. 3e), which colocalized completely (PPH4 = 1). The minor allele increased both accessibility and expression, consistent with a positive regulatory relationship and with the peak–gene SMR estimate (bSMR = 0.47, PSMR = 2.38 × 10−38). The framework also recovered the ETS2 lead eQTL (21:39094644:G:A) in a CD14Mono credible set (PIP = 0.812; Supplementary Fig. 24), matching previously prioritized functional results40 that GWAS-only fine-mapping could not resolve because of high linkage disequilibrium among candidate SNPs. By jointly modelling chromatin accessibility, gene expression and disease risk, this approach resolves the coordinated regulatory effects through which individual variants act, including at loci where high linkage disequilibrium compromises GWAS-only fine-mapping.caQTLs vary along cell-state trajectoryTo test whether genetic effects on chromatin accessibility vary with cell state, we inferred pseudotime across the CD4+ T cell compartment and tested caQTL-pseudotime interactions at single-cell resolution with CASTIE31,43 (Methods). We identified 209,189 significant interactions (false discovery rate