Genetic variants shape single-cell chromatin accessibility and disease risk

Sample collection and library processing

To build a chromatin accessibility atlas for human PBMCs, we collected blood samples from 1,042 consenting individuals from the TenK10K phase I cohort⁵, which is a combination of 952 from the Tasmanian Ophthalmic Biobank (TOB) cohort, 64 from the BioHeart (BH) cohort⁵⁵, and 26 from the Liquid Biopsy Biobank (LBIO) cohort. We obtained the 952 donors from a subset of the 982 donors of European ancestry in OneK1K¹⁰, the pilot cohort of TenK10K. We processed TOB donor blood samples in different batches than those used to process scRNA-seq samples in Cuomo et al.⁵. The BioHEART cohort was recruited between 29 March 2016 and 2 May 2023, and included participants who presented for clinically indicated cardiovascular imaging or with acute ST-elevation myocardial infarction. The LBIO was established to collect peripheral blood from participants with advanced cancer receiving intravenous therapy. Recruitment for LBIO commenced in December 2019 and is ongoing. No statistical methods were used to predetermine sample size. The study did not involve randomization or blinding.

Ethics oversight

This study was conducted in accordance with all relevant ethical regulations. It was overseen by the Human Research Ethics Committees at the Royal Children’s Hospital, Northern Sydney Local Health District, St. Vincent’s Hospital and the University of Tasmania. We obtained ethics approval for all three datasets: TOB (ethics protocol 2020/ETH02479), BioHEART study (ethics protocol 2019/ETH08376) and LBIO (2019/ETH13113). All participants provided informed consent.

Single-cell ATAC-seq

In total, we processed 952 individuals from the TOB cohort using the 10X Genomics Chromium Next GEM single-cell ATAC kit v2. This is the main dataset used for caQTL mapping. In brief, we thawed, washed, counted and pooled cryopreserved PBMC samples into 119 pools of 8 individuals, with equal numbers of live cells per donor. We isolated nuclei and transposed them before running each pool across two single-cell reactions, targeting either the same or different numbers of nuclei, ranging from 20,000 to 50,000. The specific number for each library can be found in Supplementary Table 1. We generated and sequenced 238 Illumina sequencing libraries at a read depth of approximately 30,000 reads per nucleus using S4 flow cell kits on the NovaSeq 6000.

Multiome sequencing

The BioHEART and LBIO samples were processed with the 10X Genomics Chromium Next GEM Single Cell Multiome ATAC + Gene Expression kit v1. The BioHEART multiome subcohort consisted of 64 individuals split into 4 pools of 16 individuals each, and the LBIO subcohort consisted of 26 individuals divided into 4 pools of 6–7 individuals each. As with the ATAC-only reaction, we thawed, washed, counted and pooled the cryopreserved PBMC samples in equal numbers of live cells per donor before isolating and transposing nuclei. We then ran each pooled cell suspension in a single reaction. We generated ATAC and gene expression libraries from the same pool of pre-amplified transposed DNA and cDNA. We sequenced ATAC and RNA libraries independently at approximately 30,000 reads per nucleus using S4 flow cell kits on the NovaSeq 6000.

Whole-genome sequencing and genotyping

We extracted DNA for whole-genome sequencing (WGS) from blood samples using the Qiagen DNeasy Blood and Tissue Kit, and measured DNA concentrations with Qubit. We prepared libraries using the KAPA PCR-free kit and sequenced them on an Illumina NovaSeq 6000 at 30× coverage. Cuomo et al.⁵ have provided more details on quality control and processing. We obtained genotype array data for all 952 TOB individuals from Yazar et al.¹⁰ and used them only for demultiplexing. Because 20 TOB individuals lack WGS data, we excluded them from caQTL mapping.

We subsetted common variants, including SNPs and indels (MAF ≥ 5%), from the WGS data for common-variant caQTL mapping. Variants with minor allele count < 20 or missing genotype rate > 15% were also excluded. We extracted rare variants (MAF < 5%) to perform rare-variant caQTL mapping. We used genotype array data generated with the Illumina Infinium Global Screening Array in the OneK1K study¹⁰ for demultiplexing.

Quality control and preprocessing

Read alignment (Cell Ranger ATAC)

Cell Ranger ATAC (v2.2.0) was used to align sequencing reads to the reference genome (GRCh38 v1.2.0); ‘the-force-cells’ parameter was set to 20,000, 40,000 or 50,000, depending on the targeted cell capture number for each sequencing library, resulting in 9,457,411 nuclei in total (Supplementary Table 1). Using an older version without the-force-cells function (v2.1.0), which relied on the automatic cell calling algorithm, resulted in 5,639,779 nuclei (Supplementary Table 2). We performed initial peak calling using the Cell Ranger built-in MACS2 method⁵⁶. We used Cell Ranger ARC to process the multiome sequencing libraries (RNA + ATAC) for the BioHeart and LBIO cohorts, yielding a total of 105,887 nuclei.

Demultiplexing and doublet detection

To assign nuclei to each donor and identify doublets, we used Vireo implemented with Demuxafy (v2.1.0)⁵⁷. Because we used the-force-cells parameter during Cell Ranger count, many low-quality nuclei were classified as unassigned during demultiplexing and removed (1,155,771). In the TOB cohort, we detected 3,644,297 doublets. The average doublet rate was 37.0% across all 238 libraries, compared with 32.7% if no force cell was used (Supplementary Table 3 and Supplementary Fig. 30). There are several big discrepancies in the quality control metrics between the two versions but such discrepancies are direct and expected consequences of the-force-cells setting used in the new pipeline and does not reflect any inconsistency in the underlying data. Critically, both pipelines converge on essentially the same high-quality nuclei after our standard quality control step, and we performed all downstream analyses on this quality control-filtered set. Two donors (10,003 nuclei) were not mapped to any donor genotype; thus, we retained 4,586,910 singlets from 950 donors after this step. For the multiome sequencing results from the BioHEART and LBIO cohorts, demultiplexing in the ATAC layer retained 41,248 singlets from 48 donors for BioHEART and 15,303 singlets from 12 donors for LBIO. We removed one BioHEART library and two LBIO libraries with poor sequencing quality and demultiplexing performance. In the gene expression layer, the data quality was also low, with only 11.9% of nuclei reliably assigned to a donor. Therefore, we used only the ATAC layer results for downstream analyses.

Quality controls

Nuclei that failed the following quality control thresholds were excluded from analyses: number of peak counts of more than 5,000 and less than 100,000, percentage of reads overlapping with peaks of more than 40, blacklist ratio of less than 0.01, nucleosome signal of less than 4 and TSS enrichment score of more than 4. We excluded seven sequencing libraries with a low number of nuclei after quality control (less than 300), resulting in the loss of one pool (eight donors excluded). After matching with WGS data, 20 donors could not be matched and were removed. We removed seven libraries (S0268_2, S0301_1, S0302_1, S0303_1, S0304_1, S0304_2 and S0331_1) after all quality controls. The final dataset included 3,472,552 nuclei from 922 TOB donors across 231 libraries (Supplementary Fig. 1). On average, each donor contributed 3,766 nuclei from scATAC-seq and 2,825 cells from matching scRNA-seq data.

Cell-type annotation by label transferring

To obtain fine-grained cell-type resolution for PBMCs, we adopted a two-step strategy to annotate the scATAC-seq data. We first curated a multiome bridging dataset combining eight multiome datasets (99,507 nuclei): six 10X public datasets, the BioHeart dataset and the LBIO dataset. We merged the bridging dataset into the same peak space using the ‘Reduce’ method, which merges all intersecting peaks. We then annotated it by label transfer from a CITE-seq reference (120,678 cells) provided by Hao et al.⁵⁸, and removed nuclei with a maximum prediction score of less than 0.3. In total, 95,348 nuclei (31 cell types excluding 3 non-PBMC types) were retained for the next step (Supplementary Figs. 31 and 32). After obtaining cell-type labels for the bridging multiome dataset, we built each ATAC library with anchors to the ATAC layer in the multiome bridging dataset to obtain cell-type labels. We described the number of nuclei per cell type for each library in the TOB cohort (Supplementary Table 4) and for each donor in the TOB, BioHeart and LBIO cohorts (Supplementary Tables 5 and 6). We also visualized cell-type annotations for each donor in boxplots to check extreme outliers (Supplementary Fig. 33) and in a principal component analysis (PCA) plot annotated by how many cell types fell outside the third quartile + 1.5 × interquartile range threshold (Supplementary Table 17 and Supplementary Fig. 34). The two technical replicates for the same pool also showed reasonably high concordance (Supplementary Fig. 35).

The annotation strategy was designed to maximize both annotation accuracy and caQTL discovery power given the constraints of our study design. With approximately 2,825 nuclei per donor after quality control, accurate annotation of rare cell subtypes with an abundance below approximately 0.1% is not statistically reliable, as many such types would be represented by fewer than three nuclei per donor on average. This per-donor yield is partly a consequence of the pooled library design (eight donors per pool), which places a practical ceiling on nuclei recovery. Furthermore, caQTL discovery power scales with the number of nuclei per donor per cell type (Fig. 2b), so many finer-grained subtypes lack sufficient representation to support meaningful caQTL mapping. To evaluate whether a cut-off of cell-type annotation for scATAC-seq data would materially change our results, we re-ran the caQTL analysis on chromosome 1 for four cell types spanning a range of abundances (CD14_Mono, NK, B naive and T_reg cells) across four candidate thresholds (0.3, 0.5, 0.7 and 0.9). We compared the caQTL effect estimates against those from the original, uncut analysis. caQTL slope estimates are highly concordant across thresholds (Pearson R ≥ 0.88 in all conditions; Supplementary Table 18 and Supplementary Fig. 36), indicating that our downstream results are robust to this choice. Deviation from the identity line is noticeable only at the most stringent threshold (0.9). It largely reflects severe data loss in low-abundance populations (for T_reg cells, the 0.9 cut-off leaves only approximately 4 nuclei per donor, which is not analytically viable). This analysis supports our decision to use the highest-scoring assignment without an additional cut-off. Together, these constraints motivate using 28 cell types as a level of granularity that balances annotation fidelity with statistical power.

Verifying cell-type annotation with imputed scRNA-seq data

After obtaining cell-type labels, we checked canonical gene markers for each cell type. We used the gene activity score generated in the previous step for each library. We then visualized the gene activity scores of the gene markers used in the Seurat^58,59/Azimuth reference (see Code availability) for CD14_Mono and cDC1 and confirmed that the canonical markers had enriched scores in the corresponding cell populations (Supplementary Fig. 32). We also visualized the gene activity scores of the markers used for all 26 cell types in the main scATAC-seq data (Supplementary Fig. 37).

Recalling peaks within each cell type

For each cell type, we recalled peaks from the raw fragment files using the MACS3 tool⁵⁶ integrated into SnapATAC2 (ref. ⁶⁰). The total number of peaks called in each cell type ranged from 427 (CD8⁺ proliferating) to 269,371 (CD16_Mono), with a median of 170,292 (Fig. 1e). The peaks were then merged using the Reduce function from Signac⁶¹ to define a consensus list of peaks. In total, the merged object contains 440,996 unique peaks, with minimum, median and maximum widths of 201 bp, 499 bp and 5,545 bp. We matched the peak regions against two comprehensive datasets: ENCODE (v4)⁶² with 2,348,854 peaks, and a whole-body scATAC-seq atlas⁶³ with 1,154,611 peaks, as well as the 10X 10k PBMCs ATAC v2 Chromium Controller dataset (164,487 peaks and 10,246 nuclei). After matching, 44,015 peaks (9.98%) did not overlap with any of these three datasets (Supplementary Table 7). The novel peaks are probably due to the increase in sample size, so we have more power to call peaks with low accessibility, as supported by the observation that the donor-level sparsity of novel peaks is consistently higher than that of known peaks (Supplementary Fig. 12). ChromBPNet also showed lower prediction accuracy for novel peaks than known peaks (Supplementary Fig. 11).

We annotated the peak list based on functional categories curated by overlapping the pre-annotated regions. In total, we mapped 437,561 peaks to 65 functional categories. One peak could match multiple categories, and when we calculated the proportion of each category, we used all matched pairs to compute the proportion for each category. For the ENCODE v3 annotation⁶⁴, there were five main categories: promoter-like signature (pLS), proximal enhancer-like signature (pELS), distal enhancer-like signature (dELS), DNase–H3K4me3 and CTCF only.

Dimensional reduction and visualization

We used EpiScanpy⁶⁵ to obtain a dimension reduction space for the scATAC-seq data. EpiScanpy is a toolkit for analysing single-cell epigenomic data, including scATAC-seq. We combined the quality-controlled scATAC-seq data into a single AnnData object without removing the 20 donors with unmatched WGS data (3,546,117 nuclei with 440,996 peaks), then binarized the data and selected the top 50,000 peaks based on the built-in feature variability score in EpiScanpy. We then further normalized the data such that the total count of each cell equals the median of total counts for all the cells before normalization. We then performed PCA with Scanpy⁶⁶ to reduce dimensionality. The cell clustering pattern was as expected (Supplementary Fig. 38). We integrated the scATAC-seq data with scRNA-seq in TenK10K by projecting them into a co-embedding space (Extended Data Fig. 2a). One cluster showed limited overlap between the two layers. This is due to differences in the proportion of the monocyte population between the two datasets because the blood samples were processed in different batches. A more detailed explanation is provided in supplementary note 1 of the TenK10K flagship paper⁵.

The caQTL mapping

Overall, after the nucleus-level and sample-level quality controls, a total of 3,472,552 nuclei of scATAC-seq from 922 individuals, 5,438,679 cells of scRNA-seq from 1,925 individuals (including the 922 individuals with scATAC-seq) and 95,728 nuclei of multiome data from 60 individuals (independent from the other two cohorts) were used for subsequent analyses resulting in 9 million single cells or nuclei (Extended Data Fig. 1b).

Cell-type-specific caQTL mapping

To map genetic variants to the variation of chromatin accessibility across individuals, we ran an association test between genetic variants (SNPs + indels) and chromatin accessibility levels in each of the 28 cell types. We aggregated peak counts (by summing all nuclei for one donor) into a pseudobulk matrix for each cell type. We did this separately for each repeat. In each repeat, we corrected the GC content bias using the ‘withinLaneNormalization’ function in EDAseq⁶⁷. We then removed peaks with more than 95% zero counts across donors. We converted the count data to counts per million (CPM) and then performed z-score normalization for all the peaks. We merged the normalized matrices from the two repeats by averaging the z-scores for overlapping peaks. For the unique peaks, we kept the original z-scores. Finally, we re-standardized each peak across all donors. We used the standardized matrix to generate chromatin accessibility principal components. We conducted sensitivity analyses to optimize association-mapping parameters, including normalization methods, window size and number of principal components (Supplementary Figs. 39–41). Although a smaller window size, such as 10 kb, yielded more significant associations, it would miss many associations outside the 10-kb window, and the number of tested variants overlapping with other QTL datasets would be limited. To generate consistent SNP–peak pairs between caQTL discovery and downstream colocalization analysis, we used ±1 Mb as the default window size.

The association model fits multiple covariates, including sex, age, genotype principal components and chromatin accessibility principal components. The final number of donors per cell type ranges from 120 to 922. We did not test the two rarest cell types, CD8⁺ proliferating (3 donors) and innate lymphoid cells (19 donors), because of the limited donor size. We determined the number of principal components fitted in the model using the Elbow method in the R package PCAForQTL⁶⁸. Each peak was tested against all the SNPs located within the ±1-Mb window around the peak centre, and the number of SNPs tested for each peak ranged from 293 to 30,942, with a median of 6,067. For each peak, TensorQTL⁶⁹ ran 10,000 permutations to correct P values using a β-approximation strategy. We then applied peak-level correction for each chromosome and called a caQTL if the q < 0.05, and we identified 243,225 caQTLs in total. We used graphics processing unit (GPU) acceleration to speed up the association tests, which substantially reduced computational time given the large number of peaks.

Because we measured scATAC-seq twice for each pool, we also ran the caQTL association for each repeated library separately and then compared the caQTL signals between the two runs. The significance level, quantified by −log₁₀(q value), and effect sizes are largely consistent (Pearson’s correlation > 0.90) across cell types except for those rare cell types with limited power (Supplementary Fig. 42).

Defining the cell-type specificity of caQTL effects

To quantify cell-type specificity of caQTLs, we assessed different models to explain the specificity using the approach described in Cuomo et al.⁵ (Fig. 1c, Extended Data Fig. 3 and Supplementary Fig. 5). We defined five possible scenarios that can occur for a given caQTL identified at q < 5% in one cell type (cell type A) when compared pairwise with another cell type (cell type B):

• Scenario 1 (specific caPeak): the caPeak associated with the caQTL is not sufficiently accessible (less than 1% of nuclei) in cell type B.

• Scenario 2 (specific caQTL): the caPeak is sufficiently accessible in cell type B, but no caQTL is identified.

• Scenario 3 (distinct caQTL): different caQTL variants in low linkage disequilibrium (r² < 0.5) are identified for the same caPeak in cell types A and B.

• Scenario 4 (shared caQTL, concordant effect): the same caQTL variant or variants in high linkage disequilibrium (r² ≥ 0.5) are identified for the same caPeak in both cell types, with a concordant direction of allelic effect (that is, the same directional effect of the minor allele on chromatin accessibility).

• Scenario 5 (shared caQTL, discordant effect): the same caQTL variant or variants in high linkage disequilibrium (r² ≥ 0.5) are identified for the same caPeak in both cell types, but with a discordant direction of allelic effect.

Functional annotation of caQTLs

There are 243,225 caQTLs, spanning 100,483 peaks and 158,251 unique top variants. Similar to the annotation of chromatin peaks, we overlapped the caPeaks with the five cCRE categories in ENCODE v3, as well as the top variant for each peak (Fig. 1f). The annotations of caPeaks are further divided into five subgroups based on the peak width.

Fine-mapping of causal variants for caQTLs using SuSiE

We performed fine-mapping analysis for all caPeaks using SuSiE²⁸. For each peak, we considered all variants located within a 1-Mb cis-window and calculated the corresponding variant correlation matrix using PLINK (v1.9)⁷⁰. We obtained z-scores from TensorQTL summary statistics and defined the sample size N as the number of donors containing the corresponding cell type. These summary statistics were supplied to the ‘susie_rss’ function with L = 10. The analysis yielded credible sets for each peak and marginal PIPs for individual variants.

Replication of caQTLs in published and in-house datasets

To replicate our identified caQTL signals, we collected caQTL results from two published studies^20,21. Because those published studies used different cell-type annotations, we counted a caQTL as replicated if it reached significance in the closest cell type. In other words, if at least one of our cell types contained a significant caQTL that overlapped with theirs, we counted it as being replicated. We also quantified replication of the direction of caQTL effects for Benaglio et al.²¹. As peak definitions differ across studies, we counted a replication if the peak regions overlapped between datasets, which might inflate concordance estimates. Conversely, fewer than 10% of our caQTLs replicate in previous studies, as expected given the limited statistical power of those earlier cohorts. We also conducted caQTL mapping on a combined multiome dataset comprising 60 donors from two cohorts (BioHeart and LBIO) after quality control. We used the same protocol for cell-type annotation, pseudobulking and caQTL association as in the main caQTL analysis for the TOB cohort.

Per-donor scATAC reads (n = 890 donors, because only the first repeated library was used) were realigned with STAR (v2.7.11b) in WASP mode⁷¹, retaining reads tagged vW:i:1 to correct for allelic mapping bias. We obtained heterozygous sites from OneK1K phased, imputed SNP-array genotypes (MAF ≥ 0.05, imputation R² ≥ 0.3, biallelic SNPs only). Reference and alternate snATAC counts at heterozygous sites were then quantified with GATK ASEReadCounter (v4.6.2) on per-(donor, cell type) pseudobulk binary alignment map (BAM) files. For each peak, we resolved haplotype counts for each individual at heterozygous sites based on phased genotypes and then summed them across the peak interval. To run this analysis, we developed an optimized Python re-implementation of the mixQTL ASC method²⁹ and ran it using per-peak haplotype counts within a 1-Mb cis-window flanking each peak, with default mixQTL parameters. To validate lead caQTLs, we calculated sign concordance and Spearman and Pearson correlations using TensorQTL and mixQTL effect size estimates for peak-internal lead caQTLs (that is, caQTLs whose lead variant fell within their associated peak) that were also ASC measurable. We considered lead caQTLs as ASC measurable if they passed the mixQTL default thresholds (asc_cut-off = 5, n_cut-off = 15; that is, peaks that had at least 15 donors with at least 5 SNP-intersecting reads for each haplotype within the peak). We excluded cell types with fewer than 30 such peak-internal ASC-measurable lead caQTLs from headline correlation and concordance analyses to ensure stable estimates.

Mapping rare caQTLs using SAIGE-QTL

We assessed the effects of rare variants on single-cell chromatin accessibility using the set-based test implemented in SAIGE-QTL. All SNPs and indels with MAF < 5% and missingness < 15% were included. For each peak, we tested the combined effects of variants within a cis-region defined as ±100 kb around the peak centre. We used a different window size here because it maximized the number of caQTLs in a sensitivity analysis of chromosome 1 and CD4 naive cells (Supplementary Fig. 13). Following the SAIGE-QTL framework, we performed both burden tests and sequence kernel association test (SKAT). Peak-level P values were combined using the Cauchy method, similar to the SAIGE-QTL paper³¹. Consistent with the procedure applied to common variants above, we controlled for multiple testing using the q value method and considered results significant at false discovery rate < 5%.

To rule out scenarios in which rare signals are induced by common variants in linkage disequilibrium with the rare variants tested, we further performed a conditional analysis that conditions on nearby significant common caQTLs. For each peak, we conditioned on the common variant with the lowest P value within a 100-kb cis-window centred on the peak, as identified from the TensorQTL results. Following stepwise conditional analysis, 53,195 significant rare caQTLs remained significant in the set-based test.

Regulatory effect prediction using ChromBPNet

Processing of single-cell ATAC-seq data

We first normalized the sequencing depth and aggregated the 238 scATAC-seq libraries. Then, we created pseudobulk fragment files by pooling all fragments from each predicted cell type. Next, we randomly split the fragments into two pseudo-replicates and called peaks with the ENCODE ATAC-seq pipeline (v2.2.3; https://github.com/ENCODE-DCC/atac-seq-pipeline) with default parameters. We used the overlap peak set for all downstream analyses.

ChromBPNet bias model training

To select the initial hyperparameters for the bias model, we partitioned the data into train, validation (chromosomes 8 and 20) and test (chromosomes 1, 3 and 6). Next, we trained bias models in three cell types with a range of available cells, including CD14 monocytes (n = 767,874), NK cells (n = 432,705) and B naive cells (n = 152,441), with a bias threshold factor from 0.1 to 1 at 0.1 intervals. We picked the best bias model of each cell type based on three criteria: (1) the Pearson’s R between predicted and observed non-peak counts was high; (2) the Pearson’s R between predicted and observed peak counts was low, but greater than −0.3; and (3) the model only captured bias motifs (such as TN5 or DNase-I) and did not capture any real transcription factor motif in transcription factor-MODISCO-derived motifs.

Next, we trained the ChromBPNet bias-factorized model³⁰ using the best-performing bias models from each of the three cell types. We then selected the final bias model based on two criteria: the highest Pearson’s R between predicted and observed peak counts, and the lowest Pearson’s R between predicted and observed non-peak counts, while ensuring the model’s Tn5 bias motif was corrected.

Finally, we trained the final bias model (NK cell with bias threshold factor of 1) with fivefold cross-validation. We followed the chromosome partition scheme used in the ChromBPNet method³⁰, where each chromosome is used as a test set exactly once across the folds.

ChromBPNet bias-factorized model training

We trained ChromBPNet models in 25 of 28 cell types, except CD8⁺ proliferating (4 cells), innate lymphoid cells (33 cells) and CD4⁺ proliferating (145 cells), using fivefold cross-validation with the final bias model selected in the previous step.

ChromBPNet variant effect prediction

For all top SNPs associated with a peak, we generated a 2,114-bp (ChromBPNet input length) region centred on the variant. Then, we computed cell-type-specific variant scores as the log₂ fold change in predicted total counts of the alternative versus reference allele using cell-type-specific ChromBPNet bias-factorized models. We also performed the same prediction using fine-mapped SNPs with high PIP scores in each peak, and we compared prediction accuracy across different PIP thresholds (Supplementary Fig. 10).

Transcription factor motif matching

To assess motif matches of the caQTL variants, we used the JASPAR 2026 vertebrate database⁷². For each variant, we computed contribution scores and performed motif discovery using Tangermeme (v0.5.1). We then assigned the best-matching transcription factor based on the most significant motif P value near the variant position.

Colocalized caQTL and eQTL to generate peak–gene links

To identify peak–gene links, we performed colocalization analysis using coloc³⁶ by integrating eQTL and caQTL summary statistics. We tested all significant caPeaks against genes within ±1 Mb of the peak centre. To keep a consistent variant list across colocalization analyses, we excluded indels from the summary statistics and focused only on SNP variants. The coloc results indicate the posterior probability of five different models: PPH0: no association with either trait; PPH1: association with trait 1, but not with trait 2; PPH2: association with trait 2, but not with trait 1; PPH3: association with trait 1 and trait 2, two independent SNPs; and PPH4: association with trait 1 and trait 2, one shared SNP. In our case, trait 1 is the chromatin accessibility and trait 2 is gene expression.

Although a colocalization event suggests a peak–gene association, it does not necessarily inform the directionality or effect size. To provide causal inference support for peak–gene pairs that we identified through coloc analysis, we further adopted the Mendelian randomization framework to infer causal association of the caPeaks on nearby genes using OmicsSMR⁷³. We first curated the caQTL and eQTL summary statistics into the MatrixEQTL format, and accordingly generated the BESD files as inputs for the SMR analysis. To select a significant peak for testing, we applied a nominal P value threshold < 5 × 10⁻⁸. This peak-selection threshold is much more stringent than q < 0.05, so OmicsSMR tests far fewer peak–gene pairs than coloc. For example, in CD14_Mono, coloc tested 740,024 pairs, and OmicsSMR tested 191,423 pairs. After accounting for overlap, we kept 190,494 pairs, so the OmicsSMR pairs are a subset (99.51%) of coloc pairs. When comparing the two methods, we focused only on the overlapping pairs.

Colocalization of caQTL, eQTL and GWAS

To further understand the genetic regulatory mechanisms of GWAS associations, we conducted colocalization between caQTLs and GWAS associations from 16 complex diseases using coloc, as well as 25 blood serum markers and 19 blood count traits, following the same traits tested in Margoliash et al.⁷⁴. The coloc results between eQTL and GWAS were acquired from Cuomo et al.⁵. For further comparison, we also obtained bulk eQTL summary statistics from eQTLGen⁹ (n = 31,684) and GTEx v10 whole-body tissues⁷⁵ (n = 70–816) and ran the same colocalization analysis. Similar to coloc analyses, we also ran SMR analysis between eQTL and GWAS. In the SMR analysis between eQTL and GWAS, we tested genes with a top eQTL nominal P < 5 × 10⁻⁸. We performed the same SMR analysis between caQTL and GWAS.

To assess how much GWAS signals colocalize with eQTLs and caQTLs, we used a two-step approach. First, we performed clumping analysis to select independent loci for each GWAS, using a significance threshold of P < 5 × 10⁻⁸, a linkage disequilibrium r² threshold < 0.01, and a window size ± 1 Mb. Next, we examined the raw eQTL and caQTL summary statistics, along with the outputs from coloc and SMR analyses. For each GWAS locus, we compiled a list of genes and peaks with significant QTL associations (nominal P < 5 × 10⁻⁸). A locus was considered colocalized with eQTLs or caQTLs if at least one gene or peak from the lists met the following criteria: for coloc, PPH4 > 0.8; for SMR, P_SMR < 5 × 10⁻⁸. We then aggregated results across all cell types to calculate the overall proportion of GWAS loci colocalizing (Fig. 3a,b).

Conventional coloc function assumes a single causal variant when assessing colocalization between two phenotypes. However, in practice, both phenotypes may harbour multiple causal variants. To relax this assumption, we integrated SuSiE, testing all pairs of credible sets across the two phenotypes, as implemented in the coloc.susie method⁴¹. Specifically, we performed SuSiE fine-mapping separately for each phenotype and then evaluated each pair of credible sets using conventional coloc. If any pair showed evidence of colocalization (PPH4 > 0.8), we considered the two phenotypes to colocalize.

Multi-trait fine-mapping for identifying putative causal variants

We applied mvSuSiE⁴², a variable selection method for multivariate linear regression, to identify the putative causal variants that simultaneously affect disease risk, gene expression and chromatin accessibility. Specifically, we selected genes that colocalized with GWAS (PPH4 > 0.8) as the focal genes. For each focal gene, we then selected caPeaks colocalizing with eQTL (PPH4 > 0.8) within ±1 Mb of the gene body. We then extracted summary statistics from the eQTL, caQTL mapping and GWAS and constructed a z-score matrix, with rows representing different SNPs and columns representing different phenotypes. We also set up canonical covariance matrices for the remaining phenotypes with mashr²⁷, as well as the linkage disequilibrium matrix for SNPs in the focal region with PLINK⁷⁰. We then fit an mvSuSiE model using the z-score matrix, covariance matrices and the linkage disequilibrium matrix as input.

The key output of a fitted mvSuSiE model includes a list of credible sets that are intended to capture at least one causal SNP; a PIP score quantifying the probability that a SNP is causal; and posterior estimates of SNP effects on each phenotype. We used the local false sign rate in mvSuSiE as a phenotype-level significance measurement to define the phenotypes affected by the credible sets (cut-off = 0.01). Collectively, we identified 1,085 credible sets that affect at least one layer of phenotypes (Supplementary Table 15), and we focused on credible sets that simultaneously affect disease risk, gene expression and chromatin accessibility.

Genetic regulatory interaction with cell state

Inferring pseudotime as a proxy for cell state using Palantir

To model cell-state progression, we performed pseudotime analysis using Palantir⁷⁶ on 1,176,697 cells aggregated across all pools and donors (N = 942) within the CD4 lineage (CD4 naive, CD4_TCM, CD4_TEM, CD4⁺ proliferating, CD4_CTL and T_reg cells). In brief, Palantir infers pseudotime by projecting cells into diffusion space, constructing a neighbourhood graph and estimating relative progression from a user-defined starting population, while also modelling lineage probabilities towards terminal states through a Markov chain framework. We first performed PCA using 75 principal components, followed by Harmony⁷⁷ batch correction using the sequencing pool as the batch variable. We then generated a UMAP embedding from the batch-corrected 75 principal components using a 30-nearest-neighbour graph. Next, we used these embeddings to compute 100 diffusion components, before running Palantir with k-nearest neighbour = 50 and num_waypoints = 8,000.

Testing the cell-state-dependent caQTL effect using CASTIE

Using cell-specific pseudotime estimates from Palantir, we applied CASTIE⁴³ within the SAIGE-QTL framework³¹ (v0.2.5.1), an extension of SAIGE-QTL designed to efficiently model heterogeneous genetic effects across varying cellular contexts directly from single-cell data, without pseudobulk aggregation. To reduce computational cost, we tested 493,520 cells from a subset of 386 donors in this analysis. First, we fitted a null Poisson mixed model for each gene across all cells. We included the inferred pseudotime as the dynamic (environmental) covariate, whereas sex, age, genotype principal components and chromatin-accessibility principal components were incorporated as static covariates. Next, we performed association tests for each genetic marker with MAF > 0.01 and minor allele count > 5, processing 10,000 markers per chunk. Finally, we used the ACAT test⁷⁸ to aggregate marker-level results into gene-level P values.

To visualize interaction patterns, we stratified the cell population into three pseudotime bins (early, mid and late), corresponding to pseudotime ranges of 0–0.33, 0.33–0.66 and 0.66–1.00, respectively.

Quantification of cell function scores in CD4⁺ cells

We computed cell function scores in CD4⁺ cells as cell states using scDeepID, a flexible multi-task transformer that leverages biological pathway databases as prior knowledge to train the model. First, we trained a scDeepID model with 2,390,948 CD4⁺ cells (including CD4⁺ naive, CD4_TCM, CD4_TEM, CD4_CTL, T_reg and CD4⁺ proliferating cells) from TenK10K phase I single-cell RNA-seq data. The training used expression features of 4,000 highly variable genes selected with the Scanpy function highly_variable_genes (n_top_genes = 4000, min_mean = 0.0123, max_mean = 3 and min_disp = 0.5), with the human Reactome database as a guide, filtered to the top 200 most informative pathways based on the presence of pathway genes in the highly variable gene list. We then extracted the attention-based latent space for UMAP representation and extracted cell function scores, followed by transforming each cell function score to approximate a normal distribution with zero mean and unit variance.

GRN inference using multiome data

Inferring cell-type-specific GRN with GLUE

We performed GRN inference using GLUE⁴⁷ in a cell-type-specific manner. Specifically, we aggregated cells from each major cell type across both scRNA-seq and scATAC-seq datasets and subsequently generated integrated, cell-type-specific objects. Following the GLUE framework recommendations, we adopted a standard Scanpy-based preprocessing pipeline. For scRNA-seq data, we normalized each cell to total counts across all genes, log-transformed the normalized counts and selected 6,000 highly variable genes per major cell type. We clipped each gene expression to max_value = 10, then scaled to unit variance across cells. We then applied PCA (n_pcs = 100), followed by Harmony in Scanpy, to account for technical batch effects between sequencing libraries. For scATAC-seq, we adopted the latent semantic indexing method implemented in GLUE for dimensionality reduction. We used the derived principal components to construct a nearest-neighbour graph and run UMAP in Scanpy with default parameters.

Next, we constructed the guidance graph required by GLUE by incorporating both genomic proximity and prior regulatory knowledge derived from eQTL data. For genomic distance, we connected an ATAC peak to a gene with a positive edge if it lies within a 150-kb window around the promoter of the gene. To incorporate regulatory evidence from eQTL data, we added an edge for each peak–gene pair where the peak overlaps a genomic locus associated with gene expression (eQTL). To ensure numerical stability, we incorporated self-loops into the graph. We then applied GLUE to integrate the multi-omics data using the constructed guidance graph to learn low-dimensional cell embeddings for each omics layer. We then performed cis-regulatory inference based on these embeddings. We computed the resulting GLUE regulatory score as the cosine similarity between gene and peak embeddings, with higher scores indicating stronger putative regulatory interactions between peak–gene pairs. The peak–gene pairs were annotated as ‘with support’ if at least one putative regulatory signal (such as an eQTL) colocalized with the pair. From the inferred peak–gene regulatory interactions, we utilized SCENIC⁴⁵ to build the GRN for each major cell type.

To incorporate caQTL results to enhance GRN inference, we used causal inference (SMR) results linking eQTL and caQTL to identify high-confidence peak–gene pairs with strong regulatory potential (restricted to caQTL–eQTL associations with SMR P < 5 × 10⁻⁸). We used these interactions to enrich the cis-regulatory relationships inferred by the GLUE model, which then served as input for SCENIC-based GRN construction in CD14_Mono. Integrating these curated regulatory links helps to identify transcription factor–target gene relationships that might otherwise be missed.

Comparison of inferred cis-regulatory scores between paired and unpaired multiome data

To evaluate the benefits of incorporating eQTL data and to compare differences in inferred cis-regulatory scores across datasets, we used the 10k human PBMC dataset processed with Cell Ranger ARC 2.0.0 from 10X Genomics. In parallel, we generated a randomly sampled subset of the TOB multiome cohort with a comparable size (10,000 cells) to the 10X dataset. Both datasets underwent the same preprocessing procedure, as described in the previous section. We incorporated GTEx and cell-type-specific eQTLs into GLUE model training.

Comparison of QTL-derived and multiome-derived peak–gene links

To compare peak–gene links derived from QTL colocalization with those inferred from multiome data, we performed complementary analyses in CD14 monocytes. We compared caQTL–eQTL colocalization results (PPH4) with peak–gene links inferred from multiome data using GLUE. We assessed concordance between the two approaches using Spearman correlation and further evaluated it with a binned analysis of GLUE scores across increasing levels of colocalization evidence. In addition, we repeated analyses after restricting to peak–gene pairs with colocalization evidence (PPH4 > 0.8).

Conceptually, the QTL-based framework evaluates whether a shared genetic variant influences both chromatin accessibility and gene expression, whereas GLUE infers peak–gene links from co-variation patterns across cells. These approaches capture distinct aspects of regulatory relationships and are therefore expected to yield partially overlapping but non-identical peak–gene link sets. As a technical constraint, we limited GLUE analysis to peak–gene pairs within 150 kb and to highly variable genes.

Benchmarking peak–gene link methods and GRN inference methods

To better understand the performance of the GLUE framework in peak–gene links and GRN inference, we conducted two sets of benchmarking with existing tools. First, we benchmarked the peak–gene links from GLUE⁴⁷ and pgBoost⁴⁹ on the same three validation sets (GWAS^79,80, eQTL^75,81 and CRISPR^{82,83,84,85,86,87,88,89,90}) used in pgBoost. pgBoost is a recent method that infers variant-to-gene regulatory strengths, which has demonstrated superior performance than other approaches. We modified the framework to evaluate the peak–gene links. For a given evaluation set and given method, we ranked candidate peak–gene links by linking score. At each unique linking score c assigned to any candidate link, we computed the enrichment of the set of ‘predicted’ links (linking score ≥ c) for the set of ‘true’ evaluation links as \(\frac{\frac{|{\rm{p}}{\rm{r}}{\rm{e}}{\rm{d}}{\rm{i}}{\rm{c}}{\rm{t}}{\rm{e}}{\rm{d}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{k}}{\rm{s}}\cap {\rm{t}}{\rm{r}}{\rm{u}}{\rm{e}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{k}}{\rm{s}}|}{|{\rm{t}}{\rm{r}}{\rm{u}}{\rm{e}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{k}}{\rm{s}}|}}{\frac{|{\rm{p}}{\rm{r}}{\rm{e}}{\rm{d}}{\rm{i}}{\rm{c}}{\rm{t}}{\rm{e}}{\rm{d}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{k}}{\rm{s}}|}{|{\rm{c}}{\rm{a}}{\rm{n}}{\rm{d}}{\rm{i}}{\rm{d}}{\rm{a}}{\rm{t}}{\rm{e}}\,{\rm{l}}{\rm{i}}{\rm{n}}{\rm{k}}{\rm{s}}|}}\). We defined ‘candidate links’ as all peak–gene pairs involving a gene scored by all methods and with peak-TSS distance of less than 1 Mb and more than 1 kb (for CRISPR analyses, we further restricted to peak–gene pairs experimentally tested by CRISPR). For a given evaluation set, we computed R ≤ 1, the maximum recall achieved by all methods. For each method, we then computed average enrichment across recall values (or average enrichment) as \(\frac{({\sum }_{i=1}^{{n}_{R}}{\mathrm{enrichment}}_{i}\times {(1-\mathrm{recall}}_{i}))}{{n}_{R}}\), where \({n}_{R}\) denotes the number of ‘true’ evaluation links corresponding to recall R. We obtained standard errors on average enrichment for each method, and on the difference in average enrichment between two methods, by bootstrapping genes using 1,000 bootstrap iterations.

We designated pgBoost as a reference to evaluate whether GLUE can achieve comparable results. Using the PBMC 10k multiome dataset from 10X Genomics, we extracted cells from four major cell types (CD14_Mono, CD4_TCM, CD4 naive and CD8 naive), which encompass more than 1,000 cells. We applied both GLUE (with and without TenK10K eQTL prior guidance) and pgBoost separately to the four cell types. By default, GLUE uses the top 6,000 highly variable genes for inferring peak–gene links. For benchmarking analysis, we further added the 1,949 genes from CRISPR validation datasets and resulted in 7,233 genes. This explains the difference in the number of peak–gene pairs between benchmarking using 10X multiome data and discovery analysis using TenK10K. Validation followed the same sets used in the pgBoost study, and we assessed method performance by evaluating whether predicted peak–gene pairs from each method scored higher for pairs present in the validation sets (Fig. 5a and Extended Data Fig. 9).

Second, we benchmarked the GRN inference results against CellOracle⁹¹ and SCENIC+⁴⁶, two state-of-the-art methods that can integrate single-cell multiome data to infer GRNs. We applied GLUE (with and without incorporating TenK10K eQTL and caQTL for enhancement) and the two comparison methods to the 10X Genomics PBMC dataset, deriving a list of transcription factor–target gene pairs for each approach. For quantitative evaluation, we curated three validation sets from different types of evidence: gene expression changes following transcription factor knockout⁹², chromatin immunoprecipitation followed by sequencing data^93,94,95 and functional protein–protein interactions⁹⁶. Following established benchmarking practice⁹⁷, we used EPR as the evaluation metric. Specifically, early precision is defined as the proportion of true-positive edges among the top-k-predicted edges. We computed EPR as the ratio of this early precision to the expected early precision of a random predictor, accounting for differences in network size and density. Error bars indicate 95% confidence intervals (1.96 × s.e.) computed from five independent runs using different random seeds.

Drug enrichment analysis

To assess whether single-cell caQTL data can support identifying drug targets, we performed an enrichment analysis of target–indication pairs using the Open Targets Platform⁹⁸ (OTP; v25.12), following an approach analogous to Henry et al.⁹⁹. We first used peak-to-trait SMR results to shortlist peaks with putative causal effects on selected disease traits, then integrated these with peak–gene SMR results to nominate candidate effector genes for each trait. We then intersected effector gene–trait pairs with target–indication pairs from OTP and extracted the highest clinical development phase achieved for drugs with mechanisms of action matching the relevant target–indication pair. We stratified pairs by maximum clinical phase reached (phases I, II, III, I+, II+, III+ and approved), and for each stratum constructed a 2 × 2 contingency table enumerating target–indication pairs supported or not supported by our caQTL-derived effector genes.

We note two important limitations. First, the number of target–indication pairs supported by both caQTL evidence and clinical validation (true positives) is small, so interpret the estimated odds ratios with caution. Second, and more fundamentally, the enrichment framework assumes that the absence of a caQTL-derived or eQTL-derived effector gene nomination for a given target–indication pair constitutes a negative result. However, this reflects incomplete data coverage rather than a true biological negative. Many nominally unsupported pairs may yet prove relevant as QTL mapping power improves or as additional traits are integrated. The enrichment results should therefore be interpreted as a lower-bound estimate of the concordance between QTL-nominated targets and established drug mechanisms.

Reporting summary

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