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.