Cell-type-specific eQTLs underlie the genetic architecture of complex traits

The challenge of connecting genetic effects on complex traits with their underlying gene regulatory mechanisms has motivated population-scale studies of gene expression using RNA-sequencing (RNA-seq) in bulk tissues, such as GTEx². Although these studies have identified thousands of expression quantitative trait loci (eQTLs) that contribute to variation in transcript abundance, at present, known eQTLs only explain about 10–30% of genetic effects on complex traits^5,6. This suggests that gene regulatory effects on complex traits are likely to depend on contexts that are not evident in bulk tissues, such as cell types, environments or developmental stages. Here, we sought to characterize the role of eQTL cell-type-specificity on complex traits.

Previous studies have focused on eQTLs that reach statistical significance in bulk tissues, a strategy that, due to limited power, identifies only a small fraction of suspected eQTLs. Furthermore, this strategy biases towards atypically large eQTL effects that are shared across all cell types within a tissue and are depleted in the genetic architecture of complex traits and relevant gene features³. Recently, eQTL studies using single-cell RNA-seq (scRNA-seq) have improved power to identify cell-type-specific eQTLs that are likely to contribute to gene regulatory effects on complex traits^7,8,9,10,11. However, these studies also focus on statistically significant eQTLs, leaving the overall role of eQTL cell-type specificity unclear.

We propose an alternative approach based on partitioning variation in scRNA-seq data. Rather than identifying individual eQTLs, our model unbiasedly quantifies the overall contributions of cell-type-shared and cell-type-specific eQTLs. Using this approach, we established the importance of cell-type-specific eQTLs in the genetic architecture of complex traits and partly explained why known eQTLs are depleted in functionally relevant gene features.

We set out to characterize cell-type-specific eQTLs by combining recent population-scale scRNA-seq datasets, which greatly improve cell-type resolution over bulk tissues, with statistical genetic models that unbiasedly quantify genetic effects, similar to the genomic-relatedness-based restricted maximum-likelihood (GREML) model of complex trait heritability¹².

Our approach, cell-type-informed genetic mixed-model analysis (CIGMA), is shown in Fig. 1 (Methods). For a given gene, CIGMA quantifies the variance explained by cell-type-shared eQTLs, \({\sigma }_{{\rm{g}}}^{2}\); the variance explained by eQTLs specific to cell-type c, v_c; and the overall eQTL specificity, which we define as \(\bar{v}/({\sigma }_{{\rm{g}}}^{2}+\bar{v})\), where \(\bar{v}\) averages v_c across cell types. When all v_c = 0, eQTLs are identical across cell types; otherwise, eQTLs are partly cell-type-specific, and we call the gene a cell-type-specific eGene (cs-eGene)⁹.

CIGMA has several key features for robust inference in scRNA-seq data (Methods and Supplementary Note 1). First, CIGMA accounts for cell-to-cell variation (δ) within each individual and cell type, including experimental noise. Second, CIGMA can jointly fit multiple random effects, such as cis and trans genomic regions or experimental batches. Third, our primary analyses assume that cell-type-specific eQTLs are independent across cell types for parsimony, but CIGMA can also fit a more realistic model allowing arbitrary eQTL correlations across cell types. Finally, CIGMA parameters are unbiasedly estimated with the method of moments and nonparametrically tested with jackknife.

To assess the ability of CIGMA to quantify the variance explained by cell-type-specific eQTLs, we performed simulations based on peripheral blood mononuclear cells (PBMCs) from OneK1K, which has the most individuals of available scRNA-seq datasets¹⁰. After quality control, we included 10,288 genes, 928 individuals, 7 cell types and 1,190 cells per individual on average (Methods and Supplementary Fig. 1). We applied CIGMA to real cis SNPs (within 500 kilobases (kb) of the gene body) for each gene.

We established that CIGMA is well calibrated in the absence of eQTLs by permuting genotypes across individuals (Fig. 2a). Next, we asked whether CIGMA is biased by cell-type-shared eQTLs by permuting cells across cell types. As expected, CIGMA detected the cell-type-shared eQTLs yet also provided calibrated null estimates for cell-type-specific eQTLs (Fig. 2a, Extended Data Fig. 1a and Supplementary Fig. 2). These results held across genes passing our quality control regardless of their total levels of expression or variance (Extended Data Fig. 1a). Interestingly, permuting subsets of cell types shows that eQTL specificity is reduced less by mixing cell types within the same lineage (Extended Data Fig. 1b). We found that the Wald test of CIGMA for significant cs-eGenes was calibrated when permuting individuals but slightly inflated when permuting cells, although this was negligible compared with real data signals (Supplementary Fig. 2 and Supplementary Table 1).

We next simulated data varying the number of individuals, N, and found that estimates were unbiased for all N and grew more precise as N grew (Fig. 2b and Methods). Power to detect cs-eGenes grew with N, rising from about 40% with N = 1,000 to about 100% with N = 2,000 (Supplementary Fig. 3). Power and precision also grew with the number of cells, the number of cell types and eQTL cell-type specificity (Extended Data Fig. 2, Supplementary Fig. 3 and Supplementary Table 2). We also found that CIGMA is robust to simulations that violate its assumptions on how eQTL effects are distributed across variants or cell types in theory and simulations (Extended Data Fig. 3, Supplementary Figs. 4 and 5 and Supplementary Notes 1.4 and 1.6) and to real genotype data (Supplementary Fig. 6).

We compared CIGMA with simpler methods to quantify eQTLs. First, we tested a version of CIGMA that ignores cell-to-cell variation (Methods), which resulted in about 90% deflated estimates of eQTL variance explained with realistic data, although this bias vanished in the limit of infinite cells (Fig. 2c). Second, we tested GREML¹², a common approach to quantify genetic effects on complex traits that neither distinguishes shared and specific eQTLs nor models cell-to-cell variation. With infinite cells, GREML unbiasedly estimated total eQTL variance when applied separately to each cell type (cell-type pseudobulk, CTP). However, when applied to a proxy for bulk RNA-seq that sums over cell types (overall pseudobulk, OP), GREML was biased in different directions depending on the number of cells; the mixture preserves cell-type-shared effects but partly cancels cell-type-specific effects, changing the relative strength of genetic and nongenetic variance¹³ (Methods). Third, we applied GREML to cell-type-shared and cell-type-specific components of expression by applying orthogonal tissue decomposition (OTD), which was downward biased with realistic data and remained biased with infinite cells¹⁴ (Methods). Fourth, we tested a multi-trait model, treating the expression of each cell type as a trait (BOLT-REML¹⁵), which is similar to CIGMA without cell-to-cell variation and gave similar results. Finally, we applied a GxE model to OP expression treating cell-type proportion as an ‘environment’ (GxEMM¹⁶), which also gave similar estimates on average, although they were much noisier because GxEMM uses less informative data (OP instead of CTP). GxEMM and BOLT-REML yielded biased specificity estimates due to unmodelled cell-to-cell variation (Supplementary Fig. 7). Overall, existing methods for complex traits cannot quantify cell-type-specific eQTLs in realistic scRNA-seq data because of cell-type-shared eQTLs and/or cell-to-cell variation.

We analysed 928 individuals and the seven most common cell types from the OneK1K scRNA-seq dataset¹⁰. We applied CIGMA to common cis SNPs for each of 10,288 genes, defined as SNPs within 500 kb of the gene body with minor allele frequency (MAF) greater than 5% (ref. ¹⁴). We detected cell-type-specific eQTLs for 193 genes (P < 0.05/10,288; Fig. 3b and Supplementary Fig. 8 and Supplementary Table 3), which we call statistically significant cs-eGenes. Ninety-seven of these genes do not have significant cell-type-specific residual interindividual variance, highlighting the ability of CIGMA to distinguish genetic variation. Of these 193 significant cs-eGenes from jointly testing all cell types, 64 were individually significant in two or more cell types (P < 0.05/10,288; Supplementary Fig. 9). These numbers represent conservative lower bounds as cs-eGene detection depends on power. We summarize the transcriptome-wide distribution of CIGMA estimates in Supplementary Fig. 10.

For example, one cs-eGene, CTLA4, encodes an inhibitory receptor downregulating T cell response (Fig. 3c). Cell-type-specific eQTLs in CTLA4 have previously been linked to multiple immune diseases^10,17. Another cs-eGene is SMDT1, which regulates calcium uptake in mitochondria. This gene has a known eQTL that is specific to monocytes at particular time points after pathogenic stimulation¹⁸. Our significant cs-eGenes almost perfectly overlap those detected by a fixed-effect test for cell-type-specific eQTLs in the same data, although the latter detected more genes¹¹ (Supplementary Fig. 11). We conclude that cs-eGenes can recover genes with known cell-type-specific eQTLs that are relevant to complex diseases.

We compared eQTL cell-type specificity with differential expression across cell types, a standard measure of cell-type specificity that does not consider variation across individuals. We expected these measures to partly overlap because cell-type-specific eQTLs can cause differential expression¹⁰. We found that cs-eGenes and differentially expressed genes (DEGs) have a weak but significant overlap (14 of the top 200 cs-eGenes and DEGs overlap, Spearman ρ = 0.19 between \(\bar{v}\) and variance across cell types; Supplementary Fig. 12). This confirmed that cell-type-specific eQTLs contribute to differential expression across cell types.

We next explored CIGMA’s generalized model of eQTL covariance across cell types. Although noisy at this sample size, eQTL sharing across cell types recovered known relationships (Fig. 3f). For example, the two B cell subtypes, B_IN and B_Mem, had high eQTL correlation. We found a similar pattern for residual interindividual covariances across cell types (correlation with eQTL covariance is r = 0.95; Supplementary Fig. 13), showing that cell types generally have similar responses to cis genetic and other interindividual effects. Importantly, our primary results from the simpler model of CIGMA are not biased by this more complex form of eQTL covariance (Extended Data Fig. 3 and Supplementary Note 1.4).

We jointly quantified cis and trans eQTLs, in which trans is defined to include variants on other chromosomes. Bulk RNA-seq studies have shown that cis eQTLs are more shared across tissues than trans^2,13. We found a similar pattern across cell types: cis eQTLs are mostly shared (31.1% specific, standard error (s.e.) 1.3%), whereas trans are mostly specific (59.4% specific, s.e. 4.6%; P_diff = 0.001, 1,000 permutations, one-sided test). The median gene expression heritability from trans eQTLs is 25% (s.e. 2.1%), much higher than the 3% from cis eQTLs (s.e. 0.1%; Fig. 3a and Supplementary Table 4). These results are robust to using means rather than medians (Extended Data Fig. 4) or to changing quality control thresholds (Extended Data Fig. 4 and Supplementary Fig. 14), covariate modelling assumptions (Supplementary Figs. 15–17), number of genetic PCs (Supplementary Fig. 18), pseudocount used to construct pseudobulk (Supplementary Fig. 19) or MAF filter (Supplementary Fig. 20).

Known eQTLs that are more distal from the TSS are more specific across developmental time points¹⁹ and environmental exposures²⁰. To test how eQTL specificity varies with distance to TSS, we fit CIGMA with varying cis window sizes. As expected, we found that eQTL cell-type specificity increased with window size, dropping to 30% when we reduced to 5 kb and rising to 41% when we increased to 1 Mb (Fig. 3e).

We analysed a proxy for bulk tissue RNA-seq we call OP, which sums over all cells per individual regardless of cell type. We found that cis eQTLs are 5% cell-type-specific in OP (Supplementary Fig. 21), far below the 30% specificity we observed within each cell type (Fig. 3a). This deflation of specificity is expected because bulk tissues mix cell types¹³ (Methods, Fig. 2c, Extended Data Fig. 1b, Supplementary Note 1 and Supplementary Table 5). Consequently, because trans eQTLs are more specific than cis, the trans:cis ratio is deflated in OP (83% trans) relative to each cell type (88% trans; P_diff = 0.005, 1,000 permutations, one-sided test). Overall, our results quantified how bulk tissues bias eQTLs towards cis effects that are shared across cell types.

We compared eQTL specificity with gene features that are enriched in complex traits but depleted in known eQTLs. Based on previous work, we hypothesized that genes with more specific eQTLs would have more evolutionary constraint, enhancers and gene network connections³ (Methods). Our results supported this hypothesis.

First, we tested evolutionary constraint using loss-of-function observed/expected upper bound fraction (LOEUF)²¹. We found that more-constrained genes had lower cell-type-shared and cell-type-specific cis genetic variances, consistent with bulk eQTLs^{3,4,5,13,14,22,23}. By contrast, we found that more-constrained genes had higher eQTL specificity (Fig. 4a), consistent with the hypothesis that more pleiotropic variants experience stronger selection²⁴. Second, we found that genes with higher eQTL specificity had more enhancers (Fig. 4a), consistent with observations that more complex regulatory landscapes facilitate more specific expression across tissues²⁵. Third, we examined gene network connectedness and found that shared eQTLs were less connected; as in ref. ³, this was driven by a marked enrichment in the least-connected gene decile (P = 1.9 × 10⁻⁵; Fig. 4a and Supplementary Table 6). Finally, we found that eQTL specificity decreased with average expression level and increased with gene length and differential expression across cell types (Supplementary Fig. 22). Nonetheless, adjusting for these features did not affect our conclusions (Extended Data Fig. 5), nor did changing quality control filters on CIGMA estimates (Supplementary Fig. 23), total expression (Supplementary Figs. 24 and 25), MAF (Supplementary Fig. 26), or the number of included cell types (Supplementary Fig. 27). Because bulk eQTLs are biased towards cell-type-shared effects, these results partly explain why known eQTLs are depleted in these gene features³.

We next tested these enrichments in residual interindividual variance components. We found that shared and specific residual interindividual variance shrink with constraint and grow with gene network connectivity, with no clear relationship to specificity, whereas the cell-type-specific component sharply grew with the number of enhancers (Supplementary Fig. 28). This may be explained by nongenetic and/or unmodelled trans variants (Supplementary Fig. 29), but we were underpowered to address this.

We asked what other properties distinguish genes with higher cis eQTL cell-type specificity. First, we assessed Gene Ontology (GO) terms and found that low-specificity eGenes were enriched in one term, ‘cytoplasmic translation’, whereas high-specificity eGenes were significantly enriched in 19 terms that all directly relate to immune cells²⁶ (Fig. 4c and Supplementary Table 7). These enrichments add confidence that eQTL specificity indexes cell-type-specific biology. Second, we considered drug target genes defined by OpenTargets²⁷, which were depleted in shared and specific eQTL variance without clear relationship to specificity. Third, we found the same results for genes with burden scores associated with disease risk²⁸ (Extended Data Fig. 6). Fourth, we used candidate cis-regulatory elements (cCREs)²⁹ from five immune cell types to define cCREs open in one cell type (unique) or in all five (common). As expected, the top DEGs were enriched in unique cCRE (1.2-fold, P = 0.001) and depleted in common cCRE (0.7-fold, P < 0.001; Extended Data Fig. 7 and Supplementary Table 8); cs-eGenes showed a similar but weaker pattern (respectively, 1.1-fold (P = 0.47) and 0.8-fold (P = 0.0052)), whereas cell-type-shared eGenes did not.

Having established that eQTL specificity enriches for gene features that are relevant to complex traits, we directly tested their enrichments in genetic effects on complex traits. We applied linkage disequilibrium score regression (LDSC) to test enrichment of the top 200 cs-eGenes in 10 complex traits^9,30 (Methods). We found that cs-eGenes were enriched for 4/7 blood-related traits and 0/3 less-blood-related traits (P ≤ 0.05; Fig. 4b and Supplementary Table 9). These enrichments were comparable to the enrichments for the top 200 DEGs, a standard approach to measure cell-type-specific heritability based on gene expression³⁰. We found no enrichments for complex trait heritability in the top 200 cell-type-shared eGenes, as defined by the \({\sigma }_{{\rm{g}}}^{2}\) parameter of CIGMA, or the top 200 bulk eGenes, as defined by applying GCTA to OP expression. We confirmed these results were robust to gene set size (100 or 300 genes), window size (300 kb and 700 kb, Supplementary Fig. 30 and Supplementary Table 10), and MAF filter (Supplementary Fig. 26 and Supplementary Table 10). Our primary LDSC test compares with matched random genes for robustness (Methods, Supplementary Fig. 31 and Supplementary Table 10), but we find consistent results using the standard LDSC with or without adjusting for DEGs (Supplementary Fig. 32 and Supplementary Table 10).

We next asked if CIGMA could implicate disease-causal cell types by testing disease heritability enrichment in significant cs-eGenes per each cell type (Supplementary Fig. 9). Using the same LDSC approach, we found six cell-type-trait enrichments (q < 0.1, Supplementary Table 11 and Supplementary Fig. 33). For example, we found that CD4_ET cells enrich for Crohn’s disease and ulcerative colitis (UC) heritability³¹ and that CD8_ET cells and NK cells enrich for systemic lupus erythematosus (SLE) heritability. Although CD4 cells are sharply depleted in SLE, CD8 effector T cells may play a greater part in SLE progression^32,33,34 and genetic risk^9,35.

We then used the abstract mediation model⁶ to quantify the complex trait heritability mediated by these gene sets. We found that the top 200 cs-eGenes mediate 9.0% heritability of blood-related traits, a five-fold enrichment compared with other genes (P = 0.003), whereas the top 200 DEGs, cell-type-shared eGenes and GCTA eGenes mediated only 4.9% (P = 0.004), 3.1% (P = 0.035) and 3.1% (P = 0.010), respectively (Extended Data Fig. 8 and Supplementary Table 12).

We concluded that cs-eGenes enrich for heritability in relevant complex traits, comparable to and distinct from DEGs, unlike shared eGenes.

We performed replication analyses in a second population-scale scRNA-seq dataset from California Lupus Epidemiology Study (CLUES) and Immune Variation Project (ImmVar). CLUES contains 1.2 million PBMCs from a combination of patients with SLE and healthy controls⁹. To mitigate potential confounding, we meta-analysed three subgroups: 75 patients of Asian ancestry with SLE, 70 patients of European ancestry with SLE and 70 controls of European ancestry without controls. Owing to the lower sample size, we had power only to study cis eQTLs in this dataset.

The median heritability of gene expression due to cis eQTLs was 4% in CLUES (s.e. 0.3%; Fig. 5a and Supplementary Table 13), similar to the 3% estimate in OneK1K (P_diff = 0.001, 1,000 permutations, two-sided test), but CLUES exhibited greater specificity than OneK1K (63% compared with 34%; P_diff = 0.001, 1,000 permutations, two-sided test). These results were consistent across all three CLUES subgroups (Extended Data Fig. 9). Next, as in OneK1K, we found that genes with more cell-type-specific eQTLs had greater evolutionary constraint, more enhancers and higher gene network connectivity (Fig. 5d and Supplementary Fig. 34). These results were consistent in the joint analysis across individuals in all three subgroups (Supplementary Fig. 35 and Supplementary Table 14).

We then compared eQTL estimates per gene between CLUES and OneK1K. Although CLUES had power only to identify four cs-eGenes, three of them overlap OneK1K cs-eGenes (Fig. 5b,c, Supplementary Fig. 36 and Supplementary Tables 3 and 13). More broadly, we found that OneK1K cs-eGenes were highly enriched for small cs-eGene P-values in CLUES (36/183 have P < 0.05, binomial P = 4 × 10⁻¹³), unlike the remaining genes (57/9,705 have P < 0.05, binomial P = 1). Finally, we compared eQTL estimates per gene across OneK1K, CLUES and CLUES subgroups. We found that \({\sigma }_{{\rm{g}}}^{2}\) and \(\bar{v}\) were significantly correlated across all comparisons (all P < 1 × 10⁻¹⁶; Supplementary Table 15 and Supplementary Figs. 37–39), with an average correlation of 0.27. These correlations were not significantly below 1 after accounting for estimation error (all P > 0.05/6; Methods). We conclude that eQTL specificity in PBMCs is a robust gene feature that is broadly consistent across ancestries, health states and datasets.

A leading theory to explain the limited overlap between known eQTLs and genetic effects on complex traits is that known eQTLs are primarily large effects that are shared across cell types, but complex traits are primarily driven by weak eQTLs that are cell-type specific^3,36. Our unbiased quantification of eQTLs using single-cell RNA-seq supports this hypothesis: cell-type-specific eQTLs are enriched in key gene features and complex trait heritability, unlike cell-type-shared eQTLs. This suggests that the missing gene regulation underlying complex traits, like the missing heritability of complex traits^12,37, will be primarily explained by eQTLs that are not individually detectable without massive sample sizes³⁸. We also established that eQTL cell-type specificity is a robust gene feature that is about two-fold enriched in trans compared with cis regulation and is associated with distal cis regulation, higher selective constraint, more enhancers and greater gene network connectivity.

Unlike standard eQTL studies that seek individual SNP effects, CIGMA quantifies variance explained by eQTLs in a set of SNPs. CIGMA has disadvantages: its estimates are noisier, it cannot colocalize GWAS hits^39,40,41, and it does not account for non-Gaussianity^42,43,44 or complex cell type relationships⁴⁵. But CIGMA eliminates biases from uneven eQTL detection power, which deplete known eQTLs in tissue-specificity and cell-type-specificity, enhancers, distance to TSS and fitness relevance^3,46,47. These biases from detection power are also reduced in known eQTLs that are context-specific^11,22,23, non-primary effects in a locus^48,49,50 or identified in individual cell types^9,10,51.

Our approach has limitations. First, our use of pseudobulk data ignores eQTL variation within cell types, decreasing estimates of genetic variance and specificity (Extended Data Fig. 1b). We expect that finer cell-type partitions or continuous cell states will further enrich eQTL specificity in complex trait heritability based on population genetic theory²⁴ and previous eQTL studies^{8,41,52,53,54}. Second, cell-type-specific eQTLs can depend on measurement scale⁵⁵, but this is unlikely to affect our conclusions: we confirmed our results on two scales; we controlled for differential expression across cell types; and confounding by shared eQTLs would only conservatively bias eQTL specificity. Finally, we did not study case–control differences in eQTL specificity in CLUES, which could help identify key disease genes and cell types⁵⁶.

CIGMA

We developed CIGMA to unbiasedly quantify cell-type-shared and cell-type-specific eQTLs in scRNA-seq data. CIGMA avoids bias from eQTL detection power by using a linear mixed model, similar to the GREML model of complex trait heritability¹². CIGMA models cell-type-specific pseudobulk data, which is computed by averaging across all cells in predefined cell types, one gene at a time. Mathematically, CIGMA models eQTL l in cell type c as the sum of a cell-type-shared effect (α_l) and a cell-type-specific effect (γ_lc):

Y_ic is the pseudobulk gene expression for individual i and cell type c. μ_c captures the average expression in cell type c across individuals. G_il is the genotype for individual i at eQTL l, with cell-type-shared and cell-type-specific random effects: \({\alpha }_{l}\mathop{ \sim }\limits^{\text{iid}}N(0,{\sigma }_{{\rm{g}}}^{2}/L)\) and \({\gamma }_{{lc}}\mathop{ \sim }\limits^{\text{ind}}N(0,{v}_{c}/L)\). The residual interindividual term, e, models nongenetic variation across individuals as well as unmodelled genetic variation and is also partitioned into shared and specific random effects (\({\sigma }_{e}^{2}\) and w_c)⁵⁷. Finally, ϵ_ic models variation across cells from individual i and cell type c, which may reflect experimental noise or cell subtypes and states. We use the cell-level data to precompute and subtract the cell-to-cell variation per individual and cell type, defined as the empirical variance across cells: \({\delta }_{{ic}}:= \frac{1}{{n}_{{ic}}({n}_{{ic}}-1)}{\sum }_{s=1}^{{n}_{{ic}}}{({y}_{{ics}}-{Y}_{{ic}})}^{2}\approx \mathrm{var}({{\epsilon }}_{{ic}})\), where n_ic is the number of cells for individual i in cell type c and y_ics is expression in cell s (Supplementary Note 1). In practice, we required each cell type to have more than 10 cells per individual⁵⁷. CIGMA simplifies to the additive model, GREML, if v_c, w_c and δ are 0 (ref. ¹²) and simplifies to our previous CTMM model if genetic effects are 0 (ref. ⁵⁷).

Algorithmically, CIGMA inputs genotype data and cell-type-specific pseudobulk for one gene. CIGMA outputs cell-type-shared (\({\sigma }_{{\rm{g}}}^{2}\)) and cell-type-specific (v) genetic variances, as well as shared (\({\sigma }_{e}^{2}\)) and specific (w) residual interindividual variances. We then define \(\mathrm{specificity}:= \frac{\bar{v}}{({\sigma }_{{\rm{g}}}^{2}+\bar{v})}\), where \(\bar{v}\) is the average of v_c over cell types. Approximately, specificity ≈ 1 − r_g, where r_g is the genetic correlation across cell types. We define heritability relative to total interindividual variance, \({\sigma }_{\mathrm{tot}}^{2}:= {\sigma }_{{\rm{g}}}^{2}+\bar{v}+{\sigma }_{e}^{2}+\bar{w}\), with \({h}_{\mathrm{shared}}^{2}:= \frac{{\sigma }_{{\rm{g}}}^{2}}{{\sigma }_{\mathrm{tot}}^{2}}\) and \({h}_{\mathrm{specific}}^{2}:= \frac{\bar{v}}{{\sigma }_{\mathrm{tot}}^{2}}\).

CIGMA can jointly fit multiple genotype matrices (Supplementary Note 1), such as cis and trans regions. In this case, CIGMA outputs estimates of \({\sigma }_{{\rm{g}}}^{2}\) and v for each input genotype matrix.

CIGMA can also fit a ‘Full’ model of genetic covariance across cell types by \({\gamma }_{l,}\mathop{ \sim }\limits^{\text{iid}}N(0,V)\), with \({\sigma }_{{\rm{g}}}^{2}=0\) for identification⁵⁷ (Supplementary Note 1). The simpler ‘Free’ model in equation (1) corresponds to assuming that \({V}_{{{cc}}^{{\prime} }}={\sigma }_{{\rm{g}}}^{2}\,+I\{c={c}^{{\prime} }\}{v}_{c}\) for all c and c′, that is, that cell types are independent conditional on the effect that is shared across all cell types. We show transcriptome-wide average results from the Full model in Fig. 3f, but all other main text results use the simpler Free model given the complexity and noise of scRNA-seq data at current sample sizes. Importantly, we demonstrate using theory (Supplementary Note 1.4) and simulations (Extended Data Fig. 3) that the Free model is not biased under the more realistic Full model: \({\sigma }_{{\rm{g}}}^{2}\) targets the average off-diagonal entry in V, and \({\sigma }_{{\rm{g}}}^{2}+\bar{v}\) targets the average diagonal entry.

Fitting CIGMA

We used Haseman–Elston (HE) regression to fit the parameters of CIGMA (Supplementary Note 1). HE regression is a computationally efficient and unbiased method-of-moments approach, making it suitable for aggregating inference across large-scale genomic datasets. Fixed effects were estimated using ordinary least squares.

We test the parameters of CIGMA with a Wald test using jackknife-based precision matrix estimates, as in CTMM⁵⁷. To test for cell-type-specific genetic effects, we evaluated the null hypothesis that v = 0, that is, that there are no cell-type-specific eQTLs. We used a Wald F-test with C numerator degrees of freedom and N − R denominator degrees of freedom, where C and N are the number of cell types and individuals, and R is the number of parameters in the model (including covariates). We use the same framework to test the null hypothesis of no cell-type-specific residual interindividual effects (w = 0) and to test for shared eQTLs (\({\sigma }_{{\rm{g}}}^{2}=0\)).

For completeness, we also implemented restricted maximum likelihood (REML) in R using optim and in Python using scipy, as it is more statistically efficient than HE regression (Supplementary Note 1). We incorporated Cholesky decompositions to speed up calculations, which reduced REML runtime from about 6 h to 20 min per gene with 1,000 individuals and four cell types in R. In practice, we use HE regression because it is much faster (about 30 s compared with 2 h per gene in our cis OneK1K analyses; Supplementary Table 16) and REML is biased by noisy estimates of δ (ref. ⁵⁷) (Extended Data Fig. 2h).

CIGMA partitions bulk heritability and specificity

Bulk tissues dampen cell-type-specific eQTLs by mixing cell types. The average cell-type-specific genetic variance per cell type is \(\bar{v}\), but in bulk it is \({v}_{\mathrm{bulk}}{\rm{:= }}{\sum }_{c}({\pi }_{c}^{2}+{\sigma }_{c}^{2}){v}_{c}\), where π_c and \({\sigma }_{c}^{2}\) are the mean and variance of cell-type proportion c over individuals¹³ (Supplementary Note 1). In practice, v_bulk is much lower than \(\bar{v}\), and provably so when all v_c are equal (Supplementary Note 1). As a consequence, the bulk genetic variance and specificity are deflated:

The same deflation applies to the specific residual interindividual variance, so the net effect on heritability depends on the relative specificities of genetic and residual interindividual variance.

When an additive model is fit to bulk expression, the variance explained by cell-type-specific eQTLs is further reduced (approximately \({\sum }_{c}{\pi }_{c}^{2}{v}_{c}\)) because the additive model cannot capture variation in cell-type proportion across individuals (Supplementary Note 1).

Simulations

Simulation using OneK1K data

To evaluate the performance of CIGMA, we conducted simulations on 10,228 genes from the OneK1K scRNA-seq dataset⁸. We simulated two distinct schemes: (1) permutation of genotypes across individuals, which disrupted the association between genotype and gene expression, effectively removing all eQTLs while preserving environmental effects and differences between cell types; and (2) permutation of cells across different cell types for each individual, which eliminated cell-type-specific effects—both genetic and environmental—while maintaining shared effects. We excluded genes with negative sums of genetic and residual interindividual variances (that is, \({\sigma }_{\mathrm{tot}}^{2} < 0\)), leaving 9,064 and 8,453 genes in genotype permutation and cell permutation, respectively.

Pseudobulk-level simulations

We simulated cell-type-specific pseudobulk for each individual from equation (1). Model parameters were chosen to match our estimates from the OneK1K dataset (Supplementary Note 1). We evaluated a range of scRNA-seq data parameters, including the number of individuals, cell type proportions, number of cell types, number of cells, estimation error of cell-to-cell noise, cell-type specificity and the distribution of specificity across cell types. We varied one parameter at a time, with the full list of simulated parameters provided in Supplementary Table 17. For each parameter setting, we ran 1,000 replicate simulations.

Comparison with other GREML-based methods

We fit GREML using the expectation-maximization algorithm of GCTA¹² (Supplementary Fig. 40). We applied GCTA to pseudobulk scRNA-seq data in two ways: (1) to the pseudobulk of each cell type separately (that is, columns of Y in equation (1)) and then averaging results; or (2) to OP, which averages together all cells, correcting for cell type proportion. In the OneK1K real data analysis, we corrected for the same covariates as in CIGMA, except for experimental batches, as GCTA cannot model such random effects. In Fig. 2c, GREML standard errors for heritability increase with the number of cells because the total variance decreases, making the ratio noisier.

We applied GxEMM¹⁶ to OP gene expression using its Free model and method of moments, treating cell type proportions as the ‘environment’. Like CIGMA, GxEMM partitions the genetic and nongenetic variance of gene expression into cell-type-shared and cell-type-specific components. However, GxEMM does not account for cell-to-cell variation, deflating estimates of eQTL variance, inflating estimates of nongenetic variation, and causing a complex mix of upward and downward biases in eQTL and nongenetic cell-type specificity⁵⁷ (Supplementary Fig. 7). The simple Wald test in GxEMM uses a parametric approximation to the precision matrix, which is not robust in current scRNA-seq data.

We applied the multi-trait mixed model BOLT-REML¹⁵ to the matrix of cell-type-specific pseudobulk (CTP), treating each cell type as a ‘trait’. Like Full model of CIGMA, BOLT-REML estimates general genetic and nongenetic covariance matrices across cell types, which we convert to estimates of cell-type-shared and cell-type-specific eQTL variance post hoc (Supplementary Note 1.4). BOLT-REML uses REML rather than method of moments and restricts to nonnegative variance component estimates, which adds power in complex traits but causes bias in small sample sizes⁵⁸.

We applied the orthogonal tissue decomposition framework¹⁴ to the CTP matrix to construct cell-type-shared and cell-type-specific vectors of expression. We defined shared expression as the average across cell types and then fit GREML (with GCTA) including cell-type proportions as fixed effects. We defined specific expression as the residual CTP matrix after subtracting shared expression, then applied GREML to each residualized cell type, then averaged heritability estimates over cell types.

Importantly, our results do not undermine the use of GxEMM, BOLT-REML and orthogonal tissue decomposition for their intended uses. Rather, the primary conclusion from our simulations is that standard heritability methods cannot be directly applied to current population-scale scRNA-seq datasets.

Analysis of scRNA-seq data from OneK1K

We applied CIGMA to the scRNA-seq data from PBMCs in the OneK1K cohort¹⁰. After following the quality control in ref. ¹⁰, the data span 1.27 million PBMCs from 981 individuals. Cells were classified into 14 cell types: CD4+ naive and central memory T cell (CD4_NC), CD4+ effector memory T cell (CD4_ET), CD4 + SOX4 T cell (CD4_SOX4), CD8+ naive and central memory T cell (CD8_NC), CD8+ effector memory T cell (CD8_ET), CD8 + S100B T cell (CD8_S100B), natural killer cell (NK), natural killer cell recruiting (NK_R), immature and naive B cell (B_IN), memory B cell (B_Mem), plasma cell (Plasma), classical monocyte (Mono_C) and non-classical monocyte (Mono_NC), and dendritic cell. We used SNPs with MAF > 5%, except in robustness analyses, relaxing this threshold to 1% or 0% (Supplementary Figs. 20 and 26).

For CIGMA, we conducted additional quality controls. For the one individual with technical replicates, we retained only the replicate with the largest number of cells. To approximately satisfy the Gaussian assumption in CIGMA, we required each cell type to have at least 10 cells in at least 90% of individuals, which is satisfied by seven cell types: CD4_NC, CD4_ET, CD8_ET, CD8_NC, NK, B_IN and B_Mem. We retained 10,288 autosomal genes with nonzero expression in more than 10% of individuals within each of these seven cell types. Finally, we included only the 928 individuals who each had more than 10 cells in each of these seven cell types.

To generate pseudobulk inputs to CIGMA, we normalized the total unique molecular identifier counts of each cell across all genes to 10,000 counts and then log-transformed, that is, we used log₁₀(CP10K + 1). For each gene, we computed Y_ic and δ_ic. To simplify interpretation, we scale expression such that OP has variance 1 (Supplementary Note 1); this also accounts for potential confounding due to the relationship between expression variance and evolutionary constraint⁵⁹. We adjusted for fixed effects, including sex, age, the first principal component (PC) of OP expression, and the first six genotype PCs (Supplementary Fig. 41). Cells were randomly pooled into pools in scRNA-seq. The investigators were blinded to group allocation during data collection and analysis. Experimental batches were adjusted as random effects. We fit age as a categorical covariate after binning into 5-year intervals, with additional bins for individuals younger than 25 years and older than 90 years. Genotype principal component analysis (PCA) was conducted using SNPs pruned with PLINK⁶⁰ to remove those in high linkage disequilibrium (r² > 0.2) in a sliding window of 50 SNPs and a step size of five SNPs.

We fit three different models using CIGMA. First, most of our analyses use only cis SNPs, defined as SNPs within 500 kb of the gene body. These analyses identified cs-eGenes and quantified cis heritability and specificity. Second, we extended these analyses to jointly fit cis and trans SNPs, where trans SNPs are defined as SNPs on chromosomes other than the focal gene. Third, we fit the general model of genetic covariance of CIGMA across cell types. The latter two models are noisy at current sample sizes, so we analyse only their transcriptome-wide averages. The first model is simpler and provides more robust estimates for individual genes, so we studied its results in extensive downstream analyses.

To aggregate heritability across genes, we calculated the median heritability for 9,065 genes with positive interindividual variance (\({\sigma }_{\mathrm{tot}}^{2} > 0\); Fig. 3a), as negative \({\sigma }_{\mathrm{tot}}^{2}\) renders heritability meaningless⁶¹. We confirmed these results using (1) the mean heritability rather than the median, after slightly strengthening the filter on \({\sigma }_{\mathrm{tot}}^{2}\); (2) the ratio of mean genetic variance to mean total interindividual variance; and (3) varying quality control thresholds on the standard error of heritability or total interindividual variances (Extended Data Fig. 4). To test differences in specificity between cis and trans, we performed 9,999 permutations, shuffling the heritability estimates between cis and trans for each gene.

We tested robustness of cis and trans estimates to several analytical choices. First, we varied our baseline covariate model, which treats batch as a cell-type-shared random effect and other covariates as cell-type-shared fixed effects. We found that cis, trans and environmental estimates were similar when we made either or both of these effects cell-type specific (all r > 0.99, 0.94 and 0.94, respectively; Supplementary Figs. 15–17). Second, our main analysis uses six genetic PCs to control for population stratification. We found that estimates were similar when we instead used four or eight PCs (all r > 0.99; Supplementary Fig. 18). Third, as interaction effects can depend on phenotype scale, we tested a different transformation of the underlying scRNA-seq by changing the pseudocount and found that estimates were consistent (all r > 0.95; Supplementary Fig. 19).

Analysis of scRNA-seq data from CLUES and ImmVar

We applied CIGMA to the scRNA-seq dataset from the CLUES and the ImmVar⁹. After quality control in the original study, this dataset spans 1.2 million PBMCs from 162 patients with SLE and 99 healthy controls. Cells were clustered into 11 cell types: CD14⁺ classical monocytes (cM), CD16⁺ non-classical monocytes (ncM), conventional dendritic cells (cDC), plasmacytoid dendritic cells (pDC), CD4⁺ T cells (CD4), CD8⁺ T cells (CD8), NKs, B cells (B), plasmablasts (PB), proliferating T and NKs (Prolif), and progenitor cells (Progen). For our analysis, we excluded individuals without genotype data and focused on the three largest subgroups analysed in ref. ⁹: 70 controls of European ancestry, 70 patients of European ancestry and 75 patients of Asian ancestry.

For CIGMA, we performed two types of analyses: (1) separate CIGMA analyses for each subgroup followed by meta-analysis; and (2) a mega-analysis jointly analysing all three subgroups in a single CIGMA analysis. In the first approach, we conducted quality controls independently for each subgroup using the same procedure as in OneK1K. After quality control, we retained 11,424, 10,842 and 10,786 genes in 70 European controls, 65 European patients and 74 Asian patients, respectively, involving the seven largest cell types: B, NK, CD4, CD8, cDC, cM and ncM. In CIGMA, we corrected for sex, age, cell processing batch and 10 PCs of OP expression as fixed effects and sequencing batch as a random effect. For age, similar to OneK1K, we categorized the cohort into 5-year intervals from 20 years to 70 years, with an additional group for individuals older than 70 years. Moreover, we corrected for three, four and four genotype PCs as fixed effects in European controls, European patients and Asian patients, respectively, based on elbows in the eigenvalue scree plots (Supplementary Fig. 41). Then, we conducted a meta-analysis on 10,553 genes common across all three subgroups. Cell-type-shared and cell-type-specific estimates from CIGMA were combined using inverse-variance weighting, with precision matrices estimated by jackknife resampling. For heritability and specificity, which have noisier precision estimates, we used sample-size weighting. Meta-analysed cell-type-specific genetic and residual interindividual effects were tested using Wald tests as in individual runs of CIGMA, using meta-analysed precision estimates across subgroups. In the mega-analysis, we analysed the same set of 10,553 genes across all 209 individuals. Apart from covariates used in the subgroup-based analyses, the mega-analysis adjusted for five genotype PCs, health state (case–control), cohort (CLUES–ImmVar) and ancestry.

We use Pearson correlation across genes to measure replication across cohorts. Owing to estimation error, these correlations will be below 1 even when the true underlying parameters are identical across cohorts. To model this null, we draw CIGMA pseudo-estimates from Gaussian distributions, in which the mean for each gene is its weighted average across cohorts and the standard deviation for the estimate of each cohort is given by its real data standard error. We simulate 200 datasets per gene to calculate empirical P-values.

Gene feature analysis

To investigate attributes related to cs-eGenes, we evaluated three gene-level features: LOEUF, enhancer count and connectedness in gene co-expression networks. LOEUF scores quantify the tolerance of a gene to loss-of-function mutations, serving as an approximate measure of selection strength acting on the gene. Genes with higher LOEUF scores are more tolerant to loss-of-function mutations and less conserved. We obtained LOEUF scores from the Genome Aggregation Database (gnomAD) v.2.1 (ref. ²¹). Enhancer counts reflect the regulatory complexity of a gene. We used counts from ref. ²⁵, which were derived from enhancer–gene links identified through chromatin states and the association between histone modifications and gene expression levels⁶². Gene connectedness, as defined in refs. ^3,63, was assessed by ranking genes based on their number of neighbours in co-expression networks constructed in ref. ⁶⁴. Genes with higher connectedness are likelier to have regulatory effects on more genes.

In OneK1K, we analysed 7,042 genes with positive total genetic variances (Fig. 4a). We confirmed these results using a complete set of 10,035 genes, which included genes with negative total genetic variances, and a refined subset of 6,578 genes, which excluded 464 genes with specificity s.e. exceeding 100 (Supplementary Fig. 23). In CLUES and ImmVar, we analysed 2,325 genes (Fig. 5d) after excluding genes with specificity s.e. exceeding 100 in any subgroup or total genetic variances below zero in meta-analysis. We also evaluated 3,075 genes defined solely by the latter criterion, which gave qualitatively similar results (Supplementary Fig. 34).

LD score regression

We used LD score regression (LDSC)⁶⁵ to investigate the impact of cs-eGenes on complex diseases by adapting its approach for cell-type-specific heritability enrichment³⁰. We included the top 200 genes with the most significant cell-type-specific genetic effects as defined by P-value. We defined the genomic annotation per gene as in the cis windows for CIGMA analyses (within 500 kb of the gene body). To validate our findings, we repeated the analysis using the top 100 and 300 genes and alternative window sizes of 300 kb and 700 kb.

We tested three additional gene sets for comparison: (1) shared eGenes identified by CIGMA, with the most significant cell-type-shared genetic variance; (2) additively heritable genes identified by GCTA, with the most significant OP heritability; and (3) DEGs identified by CIGMA, with the largest variance in mean expression levels (μ). The DEG set was chosen based on variance instead of significance levels because most of the genes exhibited highly significant differential gene expression across cell types⁵⁷. We chose to threshold eGenes into discrete sets, but larger datasets will enable modelling eQTL specificity as a continuous annotation.

Our analysis included seven autoimmune diseases: ulcerative colitis, rheumatoid arthritis, primary biliary cirrhosis, multiple sclerosis, SLE, Crohn’s disease and celiac disease. As negative controls, which are less relevant to immune cells, we included height, coronary artery disease and schizophrenia. GWAS summary statistics for these diseases and traits were obtained from https://alkesgroup.broadinstitute.org/sumstats_formatted.

In each analysis, LDSCs were computed using genotype data from the 1,000 Genomes Phase 3 European populations⁶⁶, restricting the analysis to SNPs in HapMap 3 and using a window size of 1 centimorgan. We removed the major histocompatibility complex (MHC) region because of its unusual LD and genetic architecture. Apart from the baseline LDSC model v.1.2 (refs. ^65,66), we conditioned on an annotation defined by all genes analysed in CIGMA so that our results are not merely biases from genes included in our study.

To control for potential confounding due to gene length and expression level, we compared LDSC P-values from each of the four tested gene sets against matched control genes to achieve empirical P-values. Control genes were selected by ranking all genes by gene length and mean expression level (average OP across individuals) and identifying, for each target gene, a matched gene that was (1) within ±500 gene ranks for both metrics and (2) located >500 kb away from all target genes. LDSC heritability enrichment analyses were then conducted on these random matched genes, replicated 999 times.

Abstract mediation model

We applied the abstract mediation model (AMM)⁶ to estimate the heritability mediated by the same gene sets (cs-eGenes, sh-eGenes, GCTA eGenes and DEGs) and complex traits as in LDSC analyses. We estimated SNP × gene rank matrix for all genes passing quality control in our OneK1K analysis, excluding the MHC region. For each SNP, we considered the closest 50 genes and binned them as suggested in ref. ⁶ (closest, 2nd closest, 3rd–5th, 6th–10th, 11th–20th, 21st–30th, 31st–40th and 41st–50th). We then estimated the fraction of heritability mediated by each bin and used these fractions to test enrichment in gene sets. Apart from individual diseases or traits, we also meta-analysed blood-related diseases (ulcerative colitis, rheumatoid arthritis, primary biliary cirrhosis, multiple sclerosis, Crohn’s disease, celiac disease and SLE) and less-blood-related diseases (height, coronary artery disease and schizophrenia).

Reporting summary

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

OneK1K single-cell gene expression and genotype data are publicly available on the Gene Expression Omnibus (GEO) under accession number GSE196830. For CLUES and ImmVar, single-cell gene expression data are available via GEO at GSE174188 and genotype data are available at dbGap under accession number phs002812.v1.p1. Summary statistics for the main CIGMA analyses in OneK1K and CLUES are provided in Supplementary Tables 3, 4 and 13. Additional publicly available data used in this study include LOEUF scores from the Genome Aggregation Database (gnomAD; https://storage.googleapis.com/gcp-public-data--gnomad/release/2.1.1/constraint/gnomad.v2.1.1.lof_metrics.by_gene.txt.bgz); enhancer counts (https://ars.els-cdn.com/content/image/1-s2.0-S0002929720300124-mmc2.xlsx); gene connectedness (https://zenodo.org/records/6618073)⁶³; GWAS summary statistics (https://alkesgroup.broadinstitute.org/sumstats_formatted); LDSC baseline annotations and 1000 Genomes Phase 3 European genotype data (https://zenodo.org/records/10515792)⁶⁶; and candidate cis-regulatory elements (cCREs; https://decoder-genetics.wustl.edu/catlasv1/humanenhancer/data/cCREs). Source data are provided with this paper.

The CIGMA Python package and all code used for analyses are available on GitHub https://github.com/Minhui-Chen/CIGMA and Zenodo (https://doi.org/10.5281/zenodo.19424343)⁶⁷.

1. Maurano, M. T. et al. Systematic localization of common disease-associated variation in regulatory DNA. Science 337, 1190–1195 (2012).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

2. GTEx Consortium. The GTEx consortium atlas of genetic regulatory effects across human tissues. Science 369, 1318–1330 (2020).

   Article  Google Scholar 

3. Mostafavi, H., Spence, J. P., Naqvi, S. & Pritchard, J. K. Systematic differences in discovery of genetic effects on gene expression and complex traits. Nat. Genet. 55, 1866–1875 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

4. Connally, N. J. et al. The missing link between genetic association and regulatory function. eLife 11, e74970 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

5. Yao, D. W., O’Connor, L. J., Price, A. L. & Gusev, A. Quantifying genetic effects on disease mediated by assayed gene expression levels. Nat. Genet. 52, 626–633 (2020).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

6. Weiner, D. J., Gazal, S., Robinson, E. B. & O’Connor, L. J. Partitioning gene-mediated disease heritability without eQTLs. Am. J. Hum. Genet. 109, 405–416 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

7. Emani, P. S. et al. Single-cell genomics and regulatory networks for 388 human brains. Science 384, eadi5199 (2024).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

8. Nathan, A. et al. Single-cell eQTL models reveal dynamic T cell state dependence of disease loci. Nature 606, 120–128 (2022).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

9. Perez, R. K. et al. Single-cell RNA-seq reveals cell type–specific molecular and genetic associations to lupus. Science 376, eabf1970 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

10. Yazar, S. et al. Single-cell eQTL mapping identifies cell type–specific genetic control of autoimmune disease. Science 376, eabf3041 (2022).

    Article  CAS  PubMed  Google Scholar 

11. Krockenberger, L. et al. FastGxC: fast and powerful context-specific eQTL mapping in bulk and single-cell data. Cell Genomics https://doi.org/10.1016/j.xgen.2026.101250 (2025).

12. Yang, J. et al. Common SNPs explain a large proportion of heritability for human height. Nat. Genet. 42, 565–569 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

13. Price, A. L. et al. Single-tissue and cross-tissue heritability of gene expression via identity-by-descent in related or unrelated individuals. PLoS Genet. 7, e1001317 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

14. Wheeler, H. E. et al. Survey of the heritability and sparse architecture of gene expression traits across human tissues. PLoS Genet. 12, e1006423 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

15. Loh, P.-R. et al. Contrasting genetic architectures of schizophrenia and other complex diseases using fast variance-components analysis. Nat. Genet. 47, 1385–1392 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

16. Dahl, A. et al. A robust method uncovers significant context-specific heritability in diverse complex traits. Am. J. Hum. Genet. 106, 71–91 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

17. Onengut-Gumuscu, S. et al. Fine mapping of type 1 diabetes susceptibility loci and evidence for colocalization of causal variants with lymphoid gene enhancers. Nat. Genet. 47, 381–386 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

18. Oelen, R. et al. Single-cell RNA-sequencing of peripheral blood mononuclear cells reveals widespread, context-specific gene expression regulation upon pathogenic exposure. Nat. Commun. 13, 3267 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

19. Strober, B. J. et al. Dynamic genetic regulation of gene expression during cellular differentiation. Science 364, 1287–1290 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

20. Lin, W. et al. Disease-associated loci share properties with response eQTLs under common environmental exposures. Preprint at bioRxiv https://doi.org/10.1101/2025.04.30.651602 (2025).

21. Karczewski, K. J. et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature 581, 434–443 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

22. GTEx Consortium. Genetic effects on gene expression across human tissues. Nature 550, 204–213 (2017).

    Article  Google Scholar 

23. Grundberg, E. et al. Mapping cis- and trans-regulatory effects across multiple tissues in twins. Nat. Genet. 44, 1084–1089 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

24. Spence, J. P. et al. Specificity, length and luck drive gene rankings in association studies. Nature 649, 918–925 (2026).

    Article  ADS  CAS  PubMed  Google Scholar 

25. Wang, X. & Goldstein, D. B. Enhancer domains predict gene pathogenicity and inform gene discovery in complex disease. Am. J. Hum. Genet. 106, 215–233 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

26. Xu, S. et al. Using clusterProfiler to characterize multiomics data. Nat. Protoc. 19, 3292–3320 (2024).

    Article  CAS  PubMed  Google Scholar 

27. Ochoa, D. et al. The next-generation Open Targets Platform: reimagined, redesigned, rebuilt. Nucleic Acids Res. 51, D1353–D1359 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

28. Karczewski, K. J. et al. Systematic single-variant and gene-based association testing of thousands of phenotypes in 394,841 UK Biobank exomes. Cell Genom. 2, 100168 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Zhang, K. et al. A single-cell atlas of chromatin accessibility in the human genome. Cell 184, 5985–6001 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

30. Finucane, H. K. et al. Heritability enrichment of specifically expressed genes identifies disease-relevant tissues and cell types. Nat. Genet. 50, 621–629 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

31. Yang, Y. et al. Investigating the shared genetic architecture between multiple sclerosis and inflammatory bowel diseases. Nat. Commun. 12, 5641 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

32. Blanco, P. et al. Increase in activated CD8+ T lymphocytes expressing perforin and granzyme B correlates with disease activity in patients with systemic lupus erythematosus. Arthritis Rheum. 52, 201–211 (2005).

    Article  CAS  PubMed  Google Scholar 

33. Buang, N. et al. Type I interferons affect the metabolic fitness of CD8+ T cells from patients with systemic lupus erythematosus. Nat. Commun. 12, 1980 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

34. McKinney, E. F., Lee, J. C., Jayne, D. R. W., Lyons, P. A. & Smith, K. G. C. T-cell exhaustion, co-stimulation and clinical outcome in autoimmunity and infection. Nature 523, 612–616 (2015).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

35. Khunsriraksakul, C. et al. Multi-ancestry and multi-trait genome-wide association meta-analyses inform clinical risk prediction for systemic lupus erythematosus. Nat. Commun. 14, 668 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

36. Farh, K. K.-H. et al. Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature 518, 337–343 (2015).

    Article  ADS  CAS  PubMed  Google Scholar 

37. Manolio, T. A. et al. Finding the missing heritability of complex diseases. Nature 461, 747–753 (2009).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

38. Yengo, L. et al. A saturated map of common genetic variants associated with human height. Nature 610, 704–712 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

39. Giambartolomei, C. et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet. 10, e1004383 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

40. Hormozdiari, F. et al. Colocalization of GWAS and eQTL signals detects target genes. Am. J. Hum. Genet. 99, 1245–1260 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

41. Mitchel, J. et al. A single-cell genetic colocalization test improves power and resolves disease-mediating cell types. Preprint at bioRxiv https://doi.org/10.1101/2025.10.10.681685 (2025).

42. Zhang, Z. E., Kim, A., Suboc, N., Mancuso, N. & Gazal, S. Efficient count-based models improve power and robustness for large-scale single-cell eQTL mapping. Preprint at medRxiv https://doi.org/10.1101/2025.01.18.25320755 (2025).

43. Zhou, W. et al. Efficient and accurate mixed model association tool for single-cell eQTL analysis. Preprint at medRxiv https://doi.org/10.1101/2024.05.15.24307317 (2024).

44. Engelmann, J. P., Palma, A., Tomczak, J. M., Theis, F. J. & Casale, F. P. Mixed models with multiple instance learning. Preprint at arxiv.org/abs/2311.02455 (2024).

45. Urbut, S. M., Wang, G., Carbonetto, P. & Stephens, M. Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions. Nat. Genet. 51, 187–195 (2019).

    Article  CAS  PubMed  Google Scholar 

46. Gamazon, E. R. et al. Using an atlas of gene regulation across 44 human tissues to inform complex disease- and trait-associated variation. Nat. Genet. 50, 956–967 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

47. Hormozdiari, F. et al. Leveraging molecular quantitative trait loci to understand the genetic architecture of diseases and complex traits. Nat. Genet. 50, 1041–1047 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

48. Dobbyn, A. et al. Landscape of conditional eQTL in dorsolateral prefrontal cortex and co-localization with schizophrenia GWAS. Am. J. Hum. Genet. 102, 1169–1184 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

49. Zeng, B. et al. Multi-ancestry eQTL meta-analysis of human brain identifies candidate causal variants for brain-related traits. Nat. Genet. 54, 161–169 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

50. Brotman, S. M. et al. Adipose tissue eQTL meta-analysis highlights the contribution of allelic heterogeneity to gene expression regulation and cardiometabolic traits. Nat. Genet. 57, 180–192 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

51. Natri, H. M. et al. Cell-type-specific and disease-associated expression quantitative trait loci in the human lung. Nat. Genet. 56, 595–604 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

52. Soskic, B. et al. Immune disease risk variants regulate gene expression dynamics during CD4+ T cell activation. Nat. Genet. 54, 817–826 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

53. Popp, J. M. et al. Cell type and dynamic state govern genetic regulation of gene expression in heterogeneous differentiating cultures. Cell Genom. 4, 100701 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

54. Qi, G. et al. Transcriptome-wide association studies at cell-state level using single-cell eQTL data. Cell Genom. 6, 101060 (2026).

    Article  CAS  PubMed  Google Scholar 

55. Ahlmann-Eltze, C. & Huber, W. Comparison of transformations for single-cell RNA-seq data. Nat. Methods 20, 665–672 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

56. Liu, X., Li, Y. I. & Pritchard, J. K. Trans effects on gene expression can drive omnigenic inheritance. Cell 177, 1022–1034 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

57. Chen, M. & Dahl, A. A robust model for cell type-specific interindividual variation in single-cell RNA sequencing data. Nat. Commun. 15, 5229 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

58. Steinsaltz, D., Dahl, A. & Wachter, K. W. Statistical properties of simple random-effects models for genetic heritability. Electron. J. Stat. 12, 321–358 (2018).

    Article  MathSciNet  PubMed  PubMed Central  Google Scholar 

59. Chen, J. et al. A quantitative framework for characterizing the evolutionary history of mammalian gene expression. Genome Res. 29, 53–63 (2019).

    Article  CAS  PubMed  Google Scholar 

60. Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4, 7 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

61. Steinsaltz, D., Dahl, A. & Wachter, K. W. On negative heritability and negative estimates of heritability. Genetics 215, 343–357 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

62. Liu, Y., Sarkar, A., Kheradpour, P., Ernst, J. & Kellis, M. Evidence of reduced recombination rate in human regulatory domains. Genome Biol. 18, 193 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

63. Mostafavi, H. Supplementary Data for ‘Systematic differences in discovery of genetic effects on gene expression and complex traits’. Zenodo. https://doi.org/10.5281/zenodo.6618073 (2023).

64. Saha, A. et al. Co-expression networks reveal the tissue-specific regulation of transcription and splicing. Genome Res. 27, 1843–1858 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

65. Finucane, H. K. et al. Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat. Genet. 47, 1228–1235 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

66. Gazal, S. S-LDSC reference files. Zenodo. https://doi.org/10.5281/zenodo.10515792 (2024).

67. Chen, M. et al. Cell type-specific eQTLs underlie the genetic architecture of complex traits. Zenodo. https://doi.org/10.5281/zenodo.19424343 (2026).

68. Wright, F. A. et al. Heritability and genomics of gene expression in peripheral blood. Nat. Genet. 46, 430–437 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

69. Gusev, A. et al. Integrative approaches for large-scale transcriptome-wide association studies. Nat. Genet. 48, 245–252 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

70. Lloyd-Jones, L. R. et al. The genetic architecture of gene expression in peripheral blood. Am. J. Hum. Genet. 100, 228–237 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

71. Liu, X. et al. Functional architectures of local and distal regulation of gene expression in multiple human tissues. Am. J. Hum. Genet. 100, 605–616 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

72. Kachuri, L. et al. Gene expression in African Americans, Puerto Ricans and Mexican Americans reveals ancestry-specific patterns of genetic architecture. Nat. Genet. 55, 952–963 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

73. Saitou, M., Dahl, A., Wang, Q. & Liu, X. Allele frequency impacts the cross-ancestry portability of gene expression prediction in lymphoblastoid cell lines. Am. J. Hum. Genet. 111, 2814–2825 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Authors and Affiliations

1. Section of Genetic Medicine, University of Chicago, Chicago, IL, USA

   Minhui Chen  (陈敏惠), Xinpei Wang  (王鑫培), Sebastian Pott, Xuanyao Liu  (刘轩尧) & Andy Dahl

2. Virginia Institute for Psychiatric and Behavioral Genetics, Virginia Commonwealth University, Richmond, VA, USA

   Minhui Chen  (陈敏惠)

3. Department of Psychiatry, Virginia Commonwealth University, Richmond, VA, USA

   Minhui Chen  (陈敏惠)

4. Bioinformatics Interdepartmental Graduate Program, University of Californaia, Los Angeles, Los Angeles, CA, USA

   Lena Krockenberger

5. Translational Genomics, Garvan Institute of Medical Research, Darlinghurst, New South Wales, Australia

   Rika Tyebally & Joseph E. Powell

6. UNSW Cellular Genomics Futures Institute, University of New South Wales, Sydney, New South Wales, Australia

   Rika Tyebally & Joseph E. Powell

7. Department of Human Genetics, University of Chicago, Chicago, IL, USA

   Jeremy J. Berg

8. Department of Human Genetics, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA, USA

   Jonathan Flint

9. Department of Pathology and Laboratory Medicine, Department of Computational Medicine, Department of Biostatistics, University of California, Los Angeles, Los Angeles, CA, USA

   Brunilda Balliu

Contributions

M.C. and A.D. developed statistical methodology, conducted analyses and wrote the paper. X.W. developed statistical methodology and conducted analyses. L.K. and R.T. assisted in data analyses. J.J.B., S.P., J.F., J.E.P., B.B. and X.L. provided suggestions and feedback to the manuscript.

Corresponding author

Correspondence to Andy Dahl.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks Shamil Sunyaev and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions