Single-nucleus transcriptome-wide association study of human brain disorders

23rd Sep 2026

Common brain disorders impose a substantial health burden, but localizing their genetic risk in the brain remains challenging¹. Although genome-wide association studies have identified numerous loci associated with neuropsychiatric and neurodegenerative disorders, many of these loci lie in non-coding regions that influence gene expression in specific cell types^2,3,4,5. Traditional bulk brain transcriptomic analyses, which often focus on European ancestry cohorts, average over cellular diversity, obscuring genetic risk-related changes in gene expression. Here we use single-nucleus gene expression profiles from the dorsolateral prefrontal cortex in the multi-ancestry PsychAD cohort to develop transcriptomic imputation models of genetically regulated expression across major brain cell types. Applying these models to neuropsychiatric and neurodegenerative disorders reveals thousands of gene–trait associations that are undetectable in bulk tissue analyses and resolves many signals to discrete neuronal, glial and immune cell populations. Cross-ancestry analyses in the Million Veteran Program confirm these associations, reveal pleiotropic effects of cell-type-specific predicted expression and demonstrate that trait-related dysregulation is conserved across ancestries, enabling mapping of causal genes and pathways. Together, these findings provide a cell-type-resolved and ancestry-aware atlas of genetically regulated expression in the human prefrontal cortex and illustrate how single-nucleus transcriptomics can sharpen gene discovery and therapeutic target prioritization for complex brain disorders.

Explore related subjects

Discover the latest articles and news in related subjects.

The genetic architecture of neuropsychiatric and neurodegenerative disorders (NPDs and NDDs, respectively) is highly polygenic, with much of their heritability attributed to common variants in non-coding regions that influence gene regulation rather than protein sequence^2,3,4,5 (Supplementary Tables 1 and 2). Genome-wide association studies (GWAS) have identified thousands of risk variants for these disorders, yet mapping these loci to effector genes and mechanisms in the human brain remains challenging. Transcriptome-wide association studies (TWAS) integrate GWAS summary statistics with predictive models of genetically regulated gene expression (GReX) to implicate genes whose expression changes may drive disease susceptibility^6,7,8,9. However, most brain TWAS reference panels are derived from homogenate bulk cortex RNA sequencing (RNA-seq) in predominantly European (EUR) ancestry cohorts^6,7,8,10, which averages over diverse neuronal and non-neuronal populations and captures only part of the genetic regulation operating in the brain. Given the extensive cellular, transcriptional^11,12 and regulatory¹³ heterogeneity of the human prefrontal cortex, approaches that treat it as a uniform tissue are poorly suited for resolving cell-type-specific contributions to NPD^2,14 and NDD⁵ risk.

The PsychAD Consortium^11,15 generated a population-scale single-nucleus RNA-seq (snRNA-seq) atlas of the dorsolateral prefrontal cortex (DLPFC), comprising over 6 million nuclei from 1,494 donors with and without major neuropsychiatric diagnoses and spanning EUR, African (AFR) and admixed American (AMR) ancestries. Leveraging this data, we constructed single-nucleus transcriptomic imputation models (snTIMs) across multiple cellular populations, enabling cell-type-resolved estimation of GReX across ancestries. We then applied snTIMs to 12 NPD and NDD GWAS to perform single-nucleus TWAS (snTWAS) and identify cell-type-specific gene–trait associations (GTAs), including previously unreported disease-linked loci. Finally, we validated and extended these findings in a large-scale phenome-wide association study (PheWAS) of approximately 600,000 participants in the Million Veteran Program (MVP). Using ancestry-matched snTIMs, we compared effect patterns across populations, characterized the pleiotropy of cell-type-specific GReX across neurological and mental health diagnoses, and refined causal gene prioritization.

snTIMs yield ancestry and cell-type-specific GReX

We developed an analytical framework to investigate how brain cell-type-specific GReX contributes to NPDs and NDDs (Extended Data Fig. 1). Using quality-controlled (see Methods) genotype and snRNA-seq data from the PsychAD Consortium^11,15, we trained 94 snTIMs across three ancestries and 32 cellular populations in the DLPFC (Extended Data Fig. 2 and Supplementary Table 3). For downstream TWAS, we restricted the analyses to snTIMs that showed robust cross-validated performance within each ancestry (see Methods). Limiting to protein-coding genes, we obtained reliable imputation models for 12,289 (74.1%), 10,375 (62.6%) and 5,543 (33.4%) of the 16,594 autosomal genes assayed in PsychAD^11,15 for EUR (N = 920), AFR (N = 321) and AMR (N = 118) ancestry, respectively (Supplementary Tables 4–9 and Supplementary Fig. 1).

The EUR single-nucleus-derived pseudobulk TIM (‘snBulk’; all nuclei pooled) and a similarly sized EUR homogenate-based TIM from the same region (‘bulk’; DLPFC EpiXcan TIM¹⁶; N = 924; 311 individuals in common with EUR snBulk) imputed a comparable number of genes (9,147 versus 8,959) with substantial but incomplete overlap (6,102 shared; 66.7–68.1%; Fig. 1a and Supplementary Note 1). Within the snTIM hierarchy, increasing the cell-type resolution from snBulk to cellular classes (‘class’) and subclasses (‘subclass’) expanded the confidently imputable set by 3,330, whereas totals remained similar between class and subclass (11,021 and 11,052, respectively; Fig. 1b, Supplementary Table S10 and Supplementary Note 2). When selecting the highest performing TIM per gene at each resolution, class and subclass snTIMs more frequently achieved the top cross-validated R² (R²_CV) than snBulk (Fig. 1b), indicating that finer cellular resolution improves model performance (Supplementary Fig. 2a). However, the presence of genes uniquely or best imputed in snBulk indicates that certain genes require larger effective sample sizes for reliable prediction.

Fig. 1: snTIM performance across cell types and ancestries.

To validate snTIM performance out of sample, we used two independent genotype–gene expression datasets: snRNA-seq from the Religious Orders Study/Memory and Aging Project (ROSMAP)¹⁷ EUR cohort and bulk RNA-seq data from fluorescence-activated cell sorting-isolated microglia (FACS-MG)^18,19. In both datasets, cross-validated and out-of-sample R² showed strong agreement, supporting the reproducible regulation of gene expression across cohorts (Supplementary Note 3).

We next compared snTIMs across ancestries. Despite differences in training sample size, shared gene–cell-type models showed strong correlations in R²_CV across ancestries (EUR–AMR: ρ = 0.58, two-sided P = 3.45 × 10^−1,158, N = 13,054; EUR–AFR: ρ = 0.54, two-sided P = 4.22 × 10^−2,711, N = 36,070; and AFR–AMR: ρ = 0.49, two-sided P = 1.29 × 10⁻⁶³³, N = 10,637; Spearman; Fig. 1c, Supplementary Table 11 and Supplementary Figs. 2b, 3 and 4), indicating that the extent to which gene expression is under genetic regulation is broadly consistent across ancestries, even when the underlying single-nucleotide polymorphisms (SNPs), which serve as gene expression predictors, differ (Supplementary Notes 4 and 5).

Together, these results show that snTIMs substantially expand the set of imputable genes compared with bulk tissue TIMs, enhance predictive performance across diverse cellular populations and capture conserved genetic regulation of brain gene expression across ancestries.

To characterize cell-type-specific GReX dysregulation in NPDs and NDDs, we applied EUR bulk and snTIMs to perform summary-level TWAS (S-TWAS) on 12 EUR GWAS (Supplementary Table 12 and external data 1 (ref. ²⁰), while controlling the false discovery rate (FDR)²¹ to account for the varying number of gene–cell-type tests at different cellular resolutions. Across traits, bulk and snBulk yielded comparable numbers of significant GTAs (1,494 versus 1,254). Among genes imputable by both TIMs, 882 were significant in snBulk and 963 in bulk, with 529 overlapping (Jaccard = 0.40). At higher levels of cellular resolution, the subclass snTIMs identified the most unique significant associations (3,003 genes across all traits), followed by the class (2,470) and snBulk (1,254) snTIMs. This progression underscores the enhanced ability of finer cellular resolution to uncover GTAs (Fig. 2a and Supplementary Fig. 5).

Fig. 2: Discovery power of S-snTWAS in individuals of EUR ancestry.

To determine whether the greater number of GTAs identified by cell-type-specific TWAS is clinically relevant, we performed gene set enrichment analysis for genes known to be associated with central nervous system-related neurological and behavioural or psychiatric symptoms in the Online Mendelian Inheritance in Man database²². Bulk and snBulk analyses revealed significant enrichment for genes linked to 5 and 4 of the 12 NPDs and NDDs, respectively, whereas class and subclass captured 11 and all 12 traits, respectively (Fig. 2b and Supplementary Table 13). This stepwise increase supports the improved biological and clinical relevance of TWAS at higher cellular resolution.

We next assessed whether higher cellular resolution preferentially identified ‘novel’ versus ‘known’ GTAs. Genes were classified as novel if the GTA was not identified by bulk S-TWAS (Supplementary Table 14; homogenate DLPFC TWAS analysis of the same GWAS) or by multi-marker analysis of genomic annotation (MAGMA) gene-based analysis^23,24 (Supplementary Table 15), which detects genes through common variant aggregation of GWAS summary statistics alone. Across all traits and snTIMs, 22.5% of significant GTAs were novel (3,711 of 16,345; Supplementary Table 16 and external data 1 (ref. ²⁰). Subclass snTIMs were more likely to detect novel than known GTAs (Fig. 2c; significant enrichment in 11 of 12 traits; Fisher’s exact test FDR-adjusted two-sided P ≤ 0.05), particularly within subclass excitatory neurons layer 6B (L6B) for bipolar disorder (BD), and subclass inhibitory neurons vasoactive intestinal peptide-expressing interneurons for anorexia nervosa (Supplementary Table 17 and Supplementary Fig. 6). These results indicate that bulk-level analyses miss a substantial fraction of cell-type-specific disease signals that are uncovered by snTWAS.

Because multiple GTAs can occur within the same locus due to linkage disequilibrium^10,25, we next performed probabilistic fine-mapping using fine-mapping of causal gene sets (FOCUS)²⁵ to identify putatively causal genes. On average, 60.8% of significant GTAs (9,936 of 16,345) were retained after fine-mapping (Fig. 2d and Supplementary Table 18), with most loci containing one gene with posterior inclusion probability (PIP) ≥ 0.5, indicating a 50% or more likelihood of being causal (Supplementary Fig. 7). Increasing cell-type resolution raised the ratio of fine-mapped genes to loci (Supplementary Note 6). This observation indicates that distinct cellular populations may contribute uniquely within shared loci, warranting further functional studies. However, lower cell-type resolution snTIMs may benefit from increased statistical power, which could contribute to improved fine-mapping precision.

S-TWAS identifies cell-type-specific disease signatures

To formally quantify the degree of cell-type specificity of GTAs within each trait, we used multivariate adaptive shrinkage (mash²⁶) to estimate posterior probabilities of cell-type-specific effects. Across all traits, 20.9% of GTAs (701 of 3,348) showed evidence of cell-type-specific effects (Supplementary Table 19), and 34.2% of these appeared in a single cell type (Fig. 3a). At the class level, class-immune exhibited the greatest degree of specificity across multiple traits (Fig. 3b, Extended Data Fig. 3a, Supplementary Table 20 and Supplementary Fig. 8). Among the top hits, CACNA1C showed highly specific effects in class-inhibitory neurons for schizophrenia (SCZ; Extended Data Fig. 3a, Supplementary Table 20 and Supplementary Fig. 8), consistent with its established functional role in SCZ pathophysiology²⁷.

Fig. 3: Cell-type specificity in snTWAS.

Existing S-TWAS methods^8,28 have yielded well-calibrated association z-scores that closely recapitulate results from individual-level TWAS results (Supplementary Table 21 and Supplementary Fig. 9a). By contrast, their effect size estimates are less reliable because they depend on reference panel quantities rather than the GWAS cohort, specifically the SNP linkage disequilibrium structure used to scale the variance of predicted expression and the allele frequencies and SNP coverage used for harmonization, and they do not propagate uncertainty in the prediction weights²⁸ (Supplementary Table 21 and Supplementary Fig. 9b,c). As a result, summary-based effect sizes are not appropriate for quantifying cell-type heterogeneity, and aggregation approaches that rely on z-scores or P values conflate true effect magnitude with sampling precision and power^29,30. To overcome this limitation, we performed individual-level snTWAS (I-snTWAS) by imputing GReX across MVP EUR. This strategy allowed us to derive NPD-specific and NDD-specific heterogeneity scores for each gene–trait pair across all imputable cell types, using the I² statistic³⁰ as a proxy (ranging from 0% to 100% for minimal to very high heterogeneity, respectively; Supplementary Table 22). In Alzheimer’s disease (AD), the top gene, BIN1, showed marked heterogeneity (I²_{cell type} = 77.36%; 95% confidence interval 78.22−91.87%; 15 imputable cell types; Fig. 3c, Supplementary Table 22 and Supplementary Fig. 10), reflecting large variation in effect size across cell types (Extended Data Fig. 3b), largely driven by cell-type-specific single-nucleus expression quantitative trait loci^12,31 (Supplementary Figs. 11 and 12). By contrast, other AD genes, including CLU, showed highly consistent effect sizes across cell types (I²_{cell type} = 0%; 95% confidence interval 0−0%; 11 imputable cell types; Fig. 3c).

We next applied a linkage disequilibrium-aware competitive pathway enrichment method³² to identify biological pathways perturbed in the summary-level snTWAS (S-snTWAS) analyses and to characterize their cell-type specificity. As cellular resolution increased, the number of significantly enriched pathways also increased (Fig. 3d, Supplementary Fig. 13 and external data 2 (ref. ²⁰). This pattern was particularly pronounced for traits such as AD and alcohol use disorder (AUD), where snTWAS signals in MG and inhibitory neurons, respectively, provided the most biologically relevant insights. For example, in AD, one of the top pathways was ‘tau protein binding’ (GO:0048156), primarily driven by dysregulation in BIN1 and APOE in class-immune (FDR-adjusted two-sided P = 1.36 × 10⁻⁶⁰) and subclass-MG (FDR-adjusted two-sided P = 2.79 × 10⁻⁶⁴; Extended Data Fig. 3c). Shared AD pathways across cell types were often driven by different genes. ‘Positive regulation of neuron death category’ (GO:1901216) was primarily driven by PICALM and APOE in class-immune (FDR-adjusted two-sided P = 8.00 × 10⁻³⁶) and subclass-MG (FDR-adjusted two-sided P = 4.49 × 10⁻³³), whereas CLU additionally contributed in other cellular populations (classes: astrocytes and oligodendrocytes; subclasses: inhibitory neuron adenosine deaminase RNA-specific B2 expressing (ADARB2)). Similarly, the ‘regulation of amyloid precursor protein catabolic process’ (GO:1902991) was driven by PICALM, APOE and SORL1 in class-immune (FDR-adjusted two-sided P = 9.73 × 10⁻⁴²) and subclass-MG (FDR-adjusted two-sided P = 7.78 × 10⁻⁴⁰), with ABCA7 further contributing in subclass-inhibitory neuron parvalbumin expressing (FDR-adjusted two-sided P = 5.83 × 10⁻³). In SCZ, which showed widespread genetically driven biological pathway dysregulation across many cellular populations (Fig. 3b), ‘divalent inorganic cation transmembrane transporter activity’ (GO:0072509) emerged as one of the strongest pathways, driven primarily by CACNA1C and GPM6A and restricted to subclass-inhibitory neuron ADARB2 (FDR-adjusted two-sided P = 1.41 × 10⁻¹³; Extended Data Fig. 3d), providing a putative cellular context for a known clinical biomarker of SCZ³³. Of note, most of these key genes were also fine-mapped in the relevant cell types (Supplementary Table 18; PIP ≥ 0.5), supporting their probable causal role.

Although pathway-level analyses revealed broad cell-type-specific perturbations, we next examined whether S-snTWAS recapitulate and extend known associations with brain pathology, using AD as a case study and focusing on MG contributions. We leveraged a hierarchical co-expression network constructed from freshly isolated human MG spanning healthy and neurodegenerative ageing¹⁹ to test MG co-expression modules for competitive enrichment of AD S-snTWAS signals (see Methods). Class-immune AD associations showed the strongest correspondence with AD pathology-linked differential expressed gene signatures¹⁹ (Extended Data Fig. 4a), a pattern not present in other cell types. Three modules (FDR-adjusted two-sided P ≤ 0.05) were significantly enriched for both AD-related differential expression and β-amyloid plaque density, and all exhibited negative enrichment (reduced predicted or observed expression in AD; Extended Data Fig. 4b and Supplementary Table 23). Functional annotation of these modules highlighted immune and myeloid activation programs, including module M824, which mapped to myeloid gene expression, and a parent–child pair (M406 → M947) pointing to TREM1 dysregulation (Extended Data Fig. 4c and Supplementary Table 24). Both M406 and M947 contained ZYX (zyxin), which our S-snTWAS predicted to be downregulated in AD (class-immune: FDR-adjusted two-sided P = 9.18 × 10⁻¹⁴, z-score = −8.29; subclass-MG: FDR-adjusted two-sided P = 1.61 × 10⁻¹³, z-score = −8.21) and fine-mapped with essentially complete posterior support (FOCUS PIP = 1.0 and 1.0; Supplementary Table 18). Of note, ZYX did not emerge in differential expression analyses (Extended Data Fig. 4d). The same modules also contained EPHA1-AS1, recently linked to AD¹⁹, which lies at the same GWAS locus as ZYX, underscoring the genetic and functional complexity of this region. Together, these findings prioritize ZYX as a putative coding effector at this locus and link its downregulation to impaired MG stress responses to β-amyloid³⁴.

Cell-type specificity reveals shared GTAs

Given the high comorbidity^35,36 and genetic correlation¹ among NPDs and NDDs, we next asked how cell-type-specific GReX dysregulation contributes to shared biology across traits. For each pair of traits, we tested the overlap of significant gene–cell-type associations using Fisher’s exact tests and found significant enrichment (FDR-adjusted two-sided P ≤ 0.05) in approximately 76% of all pairs (50 of 66; Fig. 4a and Supplementary Table 25). Comparable results were obtained when considering shared significant genes irrespective of cell type (Fig. 4a and Supplementary Table 26), indicating that these enrichments are not simply driven by cell-type redundancy. To minimize potential inflation due to overlapping samples across GWAS, we compared S-snTWAS with MVP I-snTWAS for traits with no documented sample overlap with the discovery GWAS. The direction of effects was highly consistent (90% sign concordance; 5,198 of 5,758 gene–cell-type–trait combinations; Supplementary Fig. 14). Moreover, genes shared among strongly genetically correlated disorders tend to have increasing concordance of effects at more stringent significance thresholds (Supplementary Figs. 15 and16), as illustrated in the progressive thresholding correlation analysis (PTCA; see Methods; Supplementary Fig. 17).

Fig. 4: snTWAS uncovers shared genetic effects on gene expression dysregulation across disorders.

Extending the analysis to substance use disorders further confirmed that genome-wide correlation patterns in S-snTWAS mirror those observed in GWAS (Supplementary Note 7), grouping disorders into their expected categories of NPD, NDD and substance related. Pathway-level analyses yielded a similar pattern: approximately 47.0% of trait pairs (31 of 66) shared at least one significantly enriched pathway–cell-type combinations (Supplementary Table 27 and Supplementary Fig. 18), with similar results when pathways were aggregated across all cell types (Supplementary Table 28 and Supplementary Fig. 18). Concordance analyses (Supplementary Fig. 16) and PTCA (Supplementary Fig. 17) likewise supported, on average, shared directionality of pathway dysregulation across related traits.

Cell-type resolution greatly enhanced the detection of cross-disorder convergence (Fig. 4b and Supplementary Table 29). For example, cross-disorder comparison of AD with Parkinson’s disease (PD) revealed the highest number of shared pathways in class-immune and subclass-MG (Fig. 4b), consistent with shared inflammatory mechanisms. Similarly, BD–SCZ and SCZ–major depressive disorder (MDD) comparisons showed their strongest convergence within subclass-inhibitory neuron ADARB2 and class-inhibitory neuron, respectively (Fig. 4b and Supplementary Table 29). These cell-type-specific patterns were not apparent in snBulk analyses (Fig. 4b) despite high global (for example, for SCZ–BD–MDD¹) or local (for example, for AD–PD³⁷) genetic correlation.

Together, these findings show that integrating cell-type-specific snTIMs into TWAS sharpens cross-disorder analyses, resolving convergent molecular pathways that are largely masked in bulk-level analyses.

Ancestry-specific snTWAS uncovers shared biology

Previous work has highlighted poor portability of TIMs across ancestries as a major limitation of multi-ancestry TWAS^38,39. To enable large-scale cross-ancestry interrogation of GReX dysregulation in NPDs and NDDs, we performed I-snTWAS in MVP using ancestry-matched individuals and snTIMs for eight disorders (fewer traits than in S-snTWAS due to sample size constraints and overlap between GWAS discovery samples and MVP). PTCA showed strong correlation of top-ranked GTAs across ancestries, comparable with the within-ancestry agreement between EUR S-snTWAS and EUR I-snTWAS (Fig. 5a and Supplementary Table 12). Pathway-level results from EUR and AFR I-snTWAS were likewise highly concordant (Fig. 5b). Despite the smaller effective sample size of MVP relative to the discovery GWASs (Supplementary Table 12), we identified significant GTAs in all ancestries (Supplementary Fig. 19 and external data 3 (ref. ²⁰). Given the comparatively small and heterogeneous MVP AMR sample (N = 61,073; 9.3% of the total)⁴⁰, which nevertheless showed good concordance in targeted replication⁶, we restricted in-depth I-snTWAS analyses to EUR and AFR to ensure adequate power for large-scale GReX mapping.

Fig. 5: Cross-ancestry snTWAS identify similar GTAs across ancestries.

The high I-snTWAS PTCA concordance across ancestries (Fig. 5a) indicates broadly conserved GReX dysregulation in NPDs and NDDs, yet only 8.3% of significant AFR GTAs (4 of 48) were also significant in EUR, even when aggregating across all cell types. Because replication power is limited for many of our traits in I-snTWAS, we formally evaluated the benefit of a multi-ancestry design by comparing ancestry-specific versus bi-ancestry TWAS fine-mapping. The overall fraction of significant I-snTWAS associations that fine-mapped was similar (EUR-only FOCUS: 482 of 1,545 = 31.2%; multi-ancestry FOCUS (MA-FOCUS): 493 of 1,589 = 31.0%; PIP ≥ 0.5), but the bi-ancestry fine-mapping recovered many AFR-only fine-mapped associations that were missed in the EUR-only analysis (25 versus 3; Fig. 5c,d). These patterns were consistent when collapsing to GTAs rather than trait–cell-type–gene combinations (Supplementary Fig. 20a,b), and the distribution of bi-ancestry PIPs resembled both EUR and AFR within their respective significant associations (Supplementary Fig. 20c).

We next asked whether bi-ancestry fine-mapping better prioritizes disease-relevant genes. In excitatory neurons, RHOBTB2 was not FDR significant in either EUR or AFR BD I-snTWAS (EUR z-score = −3.27, log(odds ratio (OR)) = −0.021, FDR-adjusted two-sided P = 0.21, effective N = 50,649.34; AFR z-score = −2.71, log(OR) = −0.032, FDR-adjusted two-sided P = 0.39, effective N = 14,164.83), yet the bi-ancestry inverse-variance meta-analysis was significant (z-score = −4.16, log(OR) = −0.023, FDR-adjusted two-sided P = 0.02; Fig. 5e, Supplementary Table 30 and Supplementary Fig. 21). Fine-mapping support increased from moderate within ancestries (FOCUS PIP_EUR = 0.78; PIP_AFR = 0.40) to near-certain in MA-FOCUS (PIP_MA = 0.98), prioritizing RHOBTB2 as the probable causal BD gene at this locus. The biological plausibility of RHOBTB2 is supported by functional studies showing that broad-complex, Tramtrack and Bric-à-brac (BTB)-domain variants increase excitability in human induced pluripotent stem cell-derived neurons and that altering the Drosophila sodium channel orthologue paralytic modifies RhoBTB-driven seizure-like phenotypes in vivo⁴¹. The improvement in fine-mapping parallels cross-ancestry shrinkage of credible SNP sets in multi-ancestry BD GWAS⁴², although to a lesser extent, consistent with the much smaller number of candidate genes versus SNPs per linkage disequilibrium block and the higher cross-ancestry concordance at the GReX level than at the SNP-effect level⁴³. Collectively, these results show that multi-ancestry snTWAS fine-mapping improves both cell-type-specific gene discovery and causal gene prioritization across ancestries.

NPD/NDD-associated GReX has pleiotropic effects

To investigate cell-type-specific pleiotropy of GReX, we performed a single-nucleus GReX phenome-wide association study (snGReX-PheWAS) in MVP for the top 623 significant snTWAS genes across all snTIMs (external data 4 (ref. ²⁰). We first assessed phenome-wide conservation of GReX dysregulation across ancestries and observed strong cross-ancestry agreement by PTCA (Fig. 6a) and by association effect size sign concordance (Supplementary Fig. 22). These findings replicate and extend the cross-ancestry patterns seen in our targeted I-snTWAS analyses. We then restricted the phenome scan to neurological, mental and behavioural, and sensory organ disorders within EUR ancestry (‘focused’ PheWAS) to characterize pleiotropic effects in greater detail.

In earlier I-snTWAS analyses, we showed that effect size heterogeneity across cell types, quantified by I²_{cell type}, varies by gene within a given trait (for example, AD; Fig. 3c). Using the multi-cell-type snGReX-PheWAS, we demonstrated that for a given gene (for example, CELF1), I²_{cell type} also varies across phenotypes (Fig. 6b; Supplementary Table 31 lists all significant associations for the genes included in the PheWAS). CELF1, which has no significant GTAs in snBulk, shows a wide range of I²_cell-type values across the relevant phenome where the GTAs are driven by different cellular populations (Fig. 6b and Supplementary Fig. 23). This complexity illustrates how cell-type-specific expression dysregulation shapes pleiotropic genetic effects that remain largely hidden in bulk-level analyses.

Overall, we detected significant heterogeneity among cell types in PheWAS associations for 80% of genes (355 of 444 EUR PheWAS genes imputable in more than one cell type). Thus, snGReX-PheWAS provides a systematic framework to map pleiotropic GReX effects with explicit cell-type specificity. As an example, performing a focused EUR snGReX-PheWAS of the top 20 AUD GTAs (from the S-snTWAS and I-snTWAS analyses; AUD is the trait with the largest number of GTAs in I-snTWAS (Supplementary Fig. 19); Fig. 6c and Supplementary Fig. 24) both validates the associations with AUD and reveals pleiotropy consistent with our shared pathway analyses. Across these 20 genes, the AUD snGReX association profile correlated strongly with those for tobacco use disorder (r₍₁₈₎ = 0.85, two-sided P = 1.39 × 10⁻⁵), BD (r₍₁₈₎ = 0.73, two-sided P = 2.59 × 10⁻⁴) and dementias (r₍₁₈₎ = 0.68, two-sided P = 9.09 × 10⁻⁴), in line with known AUD comorbidities and complications^44,45,46 (Fig. 6c).

Before the widespread adoption of snRNA-seq approaches, genetically regulated gene expression in human brain cell populations was examined only in low-throughput studies^5,47,48,49. More recent snRNA-seq resources have enabled a more systematic mapping of cell-type-specific expression in the human brain, but existing datasets have so far included at most several hundred EUR donors (192 (ref. ⁵⁰) and 424 (ref. ⁵¹) donors) and have not modelled GReX in other ancestries. Here we have reported a large-scale brain snTWAS that uses TIMs for EUR (N = 920), AFR (N = 321) and AMR (N = 118) ancestries to study 12 brain disorders (NPDs and NDDs), along with phenome-wide analyses in the MVP cohort (1,428 EUR, 1,057 AFR and 725 AMR phecodes). Using these snTIMs, we identified trait-associated GReX changes that are missed in homogenate tissue, resolved their cell-type specificity, characterized shared and distinct effects across disorder pairs, quantified cross-ancestry concordance and performed large-scale snGReX-PheWAS in MVP. We have provided the snTIMs, results and summary statistics as a community translational resource for target prioritization and for defining the cellular contexts in which putative risk genes act.

Class and subclass-level snTIMs imputed more genes than either pseudo-homogenate snBulk (+20% and +21%, respectively) or a similarly sized homogenate bulk model⁶ from the same region (+23% and +23%, respectively), and did so with higher prediction accuracy. In our S-snTWAS of 12 NPDs and NDDs, we found that higher cellular resolution snTIMs yield more GTAs overall, with a higher likelihood that they are novel and clinically relevant, although we note that the higher power of snBulk individually yielded more GTAs than any one other individual snTIM. Of note, S-snTWAS was FDR-adjusted across all TIMs and traits to reduce the effect of differential power. Relative to the only published brain snTWAS in MDD⁵¹, we observed 6,095 more imputable and 31 more significant gene–cell-type combinations. Our thorough investigation into cell-type specificity, using both a probabilistic framework (mash in S-snTWAS) and a frequentist heterogeneity metric (I² from MVP I-snTWAS), revealed that the most distinct signals across brain cell types arise from the class-immune population, consistent with its separate erythromyeloid origin⁵². Disease-relevant pathways were preferentially, and in some cases uniquely, detected in biologically relevant cell types.

To validate and interpret these pathway-level signals, we performed a MG-focused, module-level analysis in AD using freshly isolated human MG¹⁹. S-snTWAS module enrichments closely tracked differentially expressed gene signatures linked to AD pathology, predominantly in class-immune, and identified three negatively enriched immune modules that were also associated with amyloid plaque burden. These modules captured altered myeloid programs and TREM1 dysregulation and prioritized the ZYX–EPHA1-AS1 locus, implicating ZYX as a candidate coding effector whose predicted downregulation may impair MG stress responses to β-amyloid exposure³⁴. Together, this AD case study illustrates how cell-type-resolved TWAS bridges genetic signals to MG co-expression modules that are directly anchored to neuropathological readouts.

Beyond single-disorder analyses, snTWAS also aids in uncovering cell-type-specific cross-disorder biology. Among 12 NPDs and NDDs, we observed an increased likelihood of shared trait-associated genes and shared pathway dysregulation in 75.8% and 47.0% of trait pairwise comparisons, respectively, revealing convergent gene-level and pathway-level dysregulation that would otherwise go undetected in homogenate TWAS. For example, although snBulk TWAS identified no shared pathways between AD and PD, MG-level analyses revealed 47 shared pathways, in line with local genetic correlation evidence³⁷ and with heritability enrichment for MG annotations in both disorders^5,53. Finally, in our snGReX-PheWAS, we validated S-snTWAS findings, demonstrated that GReX cell-type heterogeneity can vary across traits for a given gene, and identified cell-type-specific associations with genetically correlated traits, comorbid conditions and trait-related complications.

In cross-ancestry analyses, snTIM performance for EUR, AFR and AMR scaled with their respective training sample sizes, as expected. Despite lower power in AFR and AMR, snTIMs in these ancestries captured an additional 882 and 335 genes that could not be reliably imputed by EUR snTIMs. On average, the degree to which the expression of a gene is under genetic control was similar across ancestries and cell types, and we observed strong cross-ancestry correlation and concordance among top-ranked GTAs in both snTWAS and snGReX-PheWAS. At the same time, non-EUR snTWAS identified ancestry-specific GTAs and improved multi-ancestry fine-mapping of risk genes. Most previous work has emphasized the poor portability of EUR-trained TIMs to other ancestries^38,39. Our work eliminates the need for mismatched ancestry application of TIMs, and we found that proper ancestry-matched study design can be highly informative for cross-ancestry NPD and NDD studies even with modest sample sizes, paralleling previous diverse ancestry blood monocyte TWASs^54,55,56.

Overall, our study extends previous work by demonstrating that cell-type-specific snTIMs are critical for identifying genetically driven, potentially actionable expression changes for translational applications. Compared with bulk approaches at the gene level, we did not observe obvious downsides to moving to snRNA-seq, aside from cost and specific scenarios where capturing extranuclear RNA is necessary. Despite the technological differences in sequencing, the sample overlap between bulk and snBulk enables their direct comparison (bulk N = 924; snBulk N = 920; 311 samples in common). Of note, we observed a similar number of imputable genes and significant S-snTWAS GTAs between bulk and snBulk, supporting this approach. For optimal model performance, future designs should consider not only increasing the number of donors but also optimizing the number of nuclei captured per specimen (Supplementary Note 4). Current snRNA-seq protocols also provide limited information on isoform abundance, even though isoform-level regulation can explain more of the heritability than gene-level expression⁵⁷. Incorporating isoform-level information is likely to further increase resolution and, in so doing, improve our understanding of the cell-type-specific architecture of brain disorders.

In summary, incorporating cell-type resolution and ancestry specificity into TWAS frameworks substantially increases gene discovery, refines causal inference and provides deeper insights into the biological underpinnings of neuropsychiatric and neurodegenerative disorders. By moving beyond bulk tissue analyses, we pave the way for more targeted therapeutic interventions that consider the intricate cellular landscape of the human brain.

Training of snTIMs

PsychAD cohort

Molecular profiling efforts of the PsychAD consortium include snRNA-seq from more than 6 million nuclei from the DLPFC of 1,494 unique donors¹⁵. As described in the Capstone paper¹¹, we utilized quadratic discriminant analysis using the 1000 Genomes Project^58,59 reference to divide the full cohort into five discrete superpopulations based on ancestry: EUR, AFR, AMR, East Asian and South Asian. Of these, 1,359 genotyped individuals (920 EUR, 321 AFR and 118 AMR)¹¹ had undergone snRNA-seq-based profiling; East Asian and South Asian were excluded from TIM building due to lack of sufficient power. Furthermore, PsychAD is composed of three ‘subcohorts’ (Mount Sinai National Institute of Health Neurobiobank, National Institute of Mental Health Intramural research program Human Brain Collection Core and Rush Alzheimer’s Disease Center (RADC)), of which RADC is predominantly AFR (Supplementary Table 41). We utilized snRNA-seq data from all 8 cell-type classes, and 23 subclasses (excluding 4 subclasses that are identical to classes). Finally, we note that we utilized release 2.5 of PsychAD.

Genotype and gene expression quality control for TIM training

Generation of an ancestry-specific common variant reference list

The merged PsychAD genotype dataset consisting of PsychAD-Mount Sinai National Institute of Health Neurobiobank SNP array, CommonMind SNP array, RADC whole-genome sequencing (WGS) and Alzheimer’s Disease sequencing project WGS samples was prepared as previously described¹⁵. Towards harmonizing SNP utilization across the training (merged PsychAD genotype dataset) and target (GWAS for S-TWAS) cohorts and trans-omics for precision medicine (TOPMed)-imputed SNP arrays for the MVP), we constructed ancestry-specific common variant reference panels. First, we queried all variants in the TOPMed-imputed^60,61,62,63 MVP’s⁶⁰ genotype release 4 and retained non-ambiguous biallelic SNPs with reference SNP cluster identification annotation included in individuals from the MVP that passed the following filters (see below for variant-level and sample-level quality control): population-specific minor allele frequency (MAF) greater than or equal to 0.01 and population-specific genotype missingness less than or equal to 0.05. Consequently, ancestry-specific SNPs were retained if they also had a matched superpopulation-specific MAF greater than or equal to 0.01 in the 1000 Genomes Project to ensure SNP overlap with GWAS used in S-snTWAS analysis^58,59. The resulting ancestry-specific common variant reference list comprises approximately 4.5, 7.5 and 5.1 million SNPs in EUR, AFR and AMR, respectively.

Variant-level filtering in the training cohort (PsychAD)

Variant-level filtering on the genotypes used for the TIM training requires that all the following variant inclusion criteria were met: SNPs exist in the ancestry-specific SNP list described above, MAF ≥ 0.01, based on ancestry-specific information from the National Center for Biotechnology Information Allele Frequency Aggregator (ALFA⁶⁴; EUR: SAMN10492695; AFR: SAMN10492698; AMR: SAMN10492700), or an in-cohort MAF ≥ 0.01 when ALFA information was not available or the variant did not pass the ALFA filter, in-sample minor allele count of 5 or greater, a Hardy–Weinberg equilibrium P value of 10⁻⁶ or greater, and not falling within areas of high linkage disequilibrium, such as the major histocompatibility complex (Supplementary Table 42). Finally, missing information for SNP predictors in existing TIMs were replaced with double the in-sample MAF (representing the in-sample average genotype), rounded to the nearest integer. Overall, we replaced about 0.050% of genotypes used for EUR snTIM training, and 0.069% each in genotypes used for AFR and AMR snTIM training.

Gene expression quality control

As reported in Lee et al.¹¹ and Zeng et al.¹², dreamlet was used to create pseudobulked gene expression by summing reads from the same individual. Expression for each cell type and each individual was computed as the log₂ counts per million with a pseudocount of 0.25 after aggregating all corresponding nuclei. A precision-weighted regression model was fit for each expressed gene in each cell type using the dreamlet package⁶⁵ using covariates for age, sex, postmortem interval, mitochondrial rate, ribosomal rate and disease status. To control for varying sequencing depth among snRNA-seq libraries, Pearson residuals (that is, residuals divided by their standard errors) were computed for each expressed gene and cell type and used in downstream analysis.

We limited genes to those within our expression annotation, Ensembl 104 (refs. ^66,67). Second, for each ‘subcohort’ within PsychAD (see ‘PsychAD cohort’), we independently scaled gene expression, performed probabilistic estimation of expression residuals (PEER)⁶⁸ to identify and adjust for hidden factors driving gene expression differences, and quantile normalized the resulting residualized gene expression. We utilized fastQTL⁶⁹ to determine the number of PEER factors (ranging from 5 to 50) that yield optimal genetic signal by identifying the point at which the number of significant expression quantitative trait loci (eQTLs; FDR (Benjamini–Hochberg method²¹)-adjusted two-sided P ≤ 0.05) is closest to 95% of the maximum number of significant eQTLs found across any assessed number of PEER factors. Finally, we combined residualized gene expression across the three subcohorts, and repeated scaling, PEER factor optimization and quantile normalization of the combined cohort to reduce the impact of the batch effect (Supplementary Fig. 30).

Training of cell-type-specific PrediXcan models

We used PrediXcan⁷ to create per-cell-type TIMs in each of three ancestries (EUR, AFR and AMR) using the quality-controlled genotypes and gene expression. Pre-processing and post-processing were done using standard scripts in R, Python 2 and Python 3. Imputable genes are considered passing R²_CV ≥ 0.01, FDR-adjusted (across cell types but within each ancestry) P_CV ≤ 0.05 and SNPs in model > 0. R²_CV and FDR-adjusted P_CV values are prediction performance R² and prediction performance P value from the ‘PredictDB’ software⁷⁰. These filtering criteria are comparable with that of other published TIMs—GTEx V8 Elastic Net models and Zeng-2024: P_CVA ≤ 0.05 (refs. ^51,71); original PrediXcan models: R²_CV > 0.01 (ref. ⁷).

To compare gene imputation models across individual TIMs or ancestries, we utilized R²_CV, a proxy of variance in gene expression explained by genetic variants serving as predictors. To compare gene imputation models across ancestries, in each pair of ancestries we selected gene–cell-type combinations confidently imputed in both ancestries, and Spearman’s correlation was performed on the resulting pairwise complete cases of confidently imputed genes.

We performed linear regression to assess the extent to which major snTIM metadata predictors (sample size and median number of nuclei contributing to the pseudobulk expression from every individual) influence snTIM performance (proxied by number of confidently imputed genes). We reported the adjusted R² as a measure of the variation explained in snTIM performance for each predictor; adjusted R² was estimated by the stats package^72,73.

We compared snTIMs against a DLPFC EpiXcan TIM¹⁶ (referred to as bulk). We compared per-gene cross-validated performance with a two-sided sign test.

Traits for S-snTWAS analysis were drawn from publicly available brain disorder GWAS. To stabilize downstream analyses, inputs were date locked in 2024. At the time of publication, the PsychAD snTIMs trained in this project will be publicly released so users can freely apply these models to newer datasets.

We began with 16 traits. Four traits showed insufficient power for S-snTWAS when analysed individually, each yielding at most one FDR-significant gene and two genome-wide significant loci in their variant association results: binge-eating disorder⁷⁴, post-traumatic stress disorder (PTSD)⁷⁵, anxiety⁷⁶ and obsessive–compulsive disorder⁷⁷. We imposed a pragmatic threshold of five genome-wide significant loci to ensure robust power for S-snTWAS analysis. Our primary analyses therefore use 12 well-powered NPD and NDD GWAS summary statistics^{75,78,79,80,81,82,83,84,85,86,87,88} (Supplementary Table 12).

For targeted comparisons of I-snTWAS versus S-snTWAS (for example, Fig. 5a), we retained anxiety⁷⁶ and PTSD⁷⁵ GWAS, as both traits had adequate power for analysis in MVP (see below). For select analyses, traits were grouped as NPDs, NDDs and substance use disorders (SUDs). These groupings are analytical and not intended to mirror formal diagnostic taxonomies. We organized traits to maximize interpretability, guided by shared genetic architecture and hypothesized dominant cell type and pathway involvement. NPDs included migraines, SCZ, BD, anorexia (anorexia nervosa), insomnia, attention-deficit/hyperactivity disorder and MDD. NDDs included AD, PD, amyotrophic lateral sclerosis and multiple sclerosis. For primary analyses, we used only AUD from the SUD group; however, for category-level comparisons, SUD representation was expanded to include tobacco use disorder and cannabis use disorder.

We utilized MungeSumstats⁸⁹ to standardize GWAS format and update reference SNP cluster identification annotations, and, where possible, we used \({n}_{\mathrm{effective}}=\,\frac{2}{\frac{1}{\mathrm{cases}}+\frac{1}{\mathrm{controls}}}\) as sample size. To ensure maximum utilization of TIM SNPs, we imputed missing variants in GWAS summary statistics based on established methods⁹⁰ (Supplementary Fig. 31). Ancestry-specific imputation linkage disequilibrium panels were built on the 1000 Genomes Project with PLINK⁹¹ and the following parameters: --ld-window-kb 1000000 --ld-window 1000 --maf 0.01 --ld-window-r2 0. We used run_imputez^90,92 using a maxWindowSize of 200. To validate imputed GWAS summary statistics, we randomly sampled 1,000 SNPs per chromosome common to the GWAS and linkage disequilibrium reference panel, and correlated real and imputed z-scores (Supplementary Fig. 31). Only SNPs with a GWAS imputation R² ≥ 0.7 were retained. Finally, we performed S-TWAS for our TIMs (external data 1 (ref. ²⁰) using summary-PrediXcan (S-PrediXcan)²⁸ and Python 2. Of note, we filtered out genes located within the major histocompatibility complex, due to high linkage disequilibrium. Significant gene–cell-type combinations were determined using the FDR-adjusted two-sided P value²¹ cut-off of 0.05. The FDR-adjusted two-sided P value was calculated across all gene–cell-type combinations for each trait. Finally, we note that we did not perform GWAS imputation for tobacco use disorder and cannabis use disorder.

Enrichment for clinically relevant genes

To validate associations from the S-TWAS, we performed gene set enrichment analysis for genes with entries in the ‘neurologicCentralNervousSystem’ and ‘neurologicBehavioralPsychiatricManifestations’ columns of the Clinical Synopsis tables of the Online Mendelian Inheritance in Man database²². For our analysis, we only considered gene-specific Online Mendelian Inheritance in Man entries and not entries corresponding to large loci (for example, regions spanning multiple genes), resulting in a final set of 1,949 genes. Query TWAS genes (nominal two-sided P ≤ 0.01; only protein-coding genes) were tested for enrichment against this gene set using a one-sided Fisher’s exact test followed by FDR adjustment²¹.

Enrichment of novel genes

We qualified a GTA as novel if the association was not present in bulk DLPFC S-TWAS (FDR-adjusted two-sided P ≤ 0.05 (ref. ²¹); FDR adjustment performed across all traits to mimic correction in S-snTWAS analysis) or MAGMA analysis (FDR-adjusted two-sided P ≤ 0.05; FDR adjustment performed across all traits to mimic correction in S-snTWAS analysis) performed using the MAGMA²⁴ SNP2GENE function on the Functional Mapping and Annotation of GWAS (FUMA) web portal^23,24. GWAS were prepared as above and then lifted over to GRCh37 using liftover within MungeSumstats^89,93 as required by FUMA^89,93. FUMA default parameters were used, except for the following: ‘genetype’ was set to ‘all’; the major histocompatibility complex region was custom defined to ‘28477797–33448354’ (to match the definition used in other analyses, accounting for genome build differences); window was set to 50 kb (upstream and downstream). The union () of FDR-significant S-TWAS and MAGMA hits was used to define our known list of GATs. Significant GTAs from the S-snTWAS not within this list were defined as novel. To test whether novel GTAs were enriched in snTIMs, we separated significant GTAs into two categories: significant in snBulk, and significant only in snTWAS (either class or subclass level; FDR-significant across all cellular populations). We used a two-sided Fisher’s exact test to obtain ORs and Bonferroni-corrected two-sided P values for the enrichment of novel GTAs in cell-type-specific TIMs.

To address potential horizontal pleiotropic effects and account for linkage disequilibrium among SNPs utilized in snTIMs, we utilized FOCUS²⁵. FOCUS models the marginal TWAS z-scores as a multivariate Gaussian distribution given the estimated eQTL effect size and the SNP correlation and utilizes a Bayesian approach to calculate the marginal PIP for each gene, indicating its likelihood of being causal in a specific TWAS risk region. Owing to internal GWAS imputation in FOCUS, we utilized post-munging GWAS summary statistics from our pipeline before missing SNP imputation. PrediXcan-based snTIMs were converted to FOCUS-compatible format after retaining only imputable genes. FOCUS was applied both individually to each snTIM to assess effects within each cell type, and jointly (multi-cell type) for further prioritization. Unless otherwise specified, we considered TWAS-significant associations with PIP ≥ 0.5 to be fine-mapped. To quantify the proportion of S-snTWAS-significant GTAs that were fine-mapped, we summed, for each trait–snTIM combination, the number of significant GTAs with fine-mapping support and divided by the total number of significant GTAs.

We evaluated whether fine-mapped genes increase in class and subclass aggregates by assessing the total number of unique fine-mapped genes and loci in each snTIM versus the respective aggregate. Because our primary analyses utilized FOCUS applied to individual snTIMs, we repeated the analysis using multi-cell-type fine-mapping. We also applied the Wilcoxon rank-sum test to compare genes per locus in aggregates versus individual snTIMs.

To perform bi-ancestry fine-mapping for EUR and AFR I-snTWAS in MVP, we used MA-FOCUS⁹⁴. We utilized TOPMed-imputed⁶¹ genotypes in MVP so that no additional imputation was needed for SNP predictors in the snTIMs; individual-level SNP missingness was handled as above. To visualize the top AFR fine-mapped associations and their concordance with EUR fine-mapped associations and MA-FOCUS, we plotted the scaled I-snTWAS z-scores within each cell type and ancestry (for meta-analysis, we performed inverse-variance weighted meta-analysis using I-snTWAS effect size and standard error and scaled the resulting z-score). We restricted the analysis to class-level snTIMs and prioritized the top ten associations by ranking each gene by the maximum scaled AFR I-snTWAS z-score × AFR fine-mapping PIP across all considered cell-types.

We implemented mash²⁶ using the mashr package in R. Owing to differences in the imputability of genes across snTIMs at the subclass level, GReX matrices are sparse; thus, we limited analysis to class-level snTIMs (Supplementary Fig. 32). Moreover, we excluded EUR classes Endo and mural, which confidently imputed a much smaller number of genes in comparison to other class-level snTIMs. Then, we performed mash on z-scores set as effect size and standard errors set to 1; we did not use the reported association effect sizes due to limitations in accurate effect size estimation in S-PrediXcan²⁸ and other TWAS methods²⁸. We filled missing values (where we do not have a model for a given gene in a given cell type), with a z-score of 0 and an arbitrarily high standard error (1 × 10¹⁰). We ran mash using data-driven covariances as recommended by the mashr authors. The mash analysis was used to obtain the probabilities of an effect (GTA) existing in a given cell type (for example, Extended Data Fig. 3a) and to assess cell-type specificity in S-snTWAS (Fig. 3a,b). For the latter, we first limited the analysis to all S-snTWAS-significant GTAs within the six cell types utilized in the mash analysis above, and further restricted GTAs to the ones having at least one but less than six (all) cell types with a significant association (local false sign rate ≤ 0.05; due to Bayesian statistics, not all S-snTWAS significant GTAs are necessarily significant in mash results).

I-snTWAS association analyses were performed as described above. For comparison, S-snTWAS association analyses were conducted using S-PrediXcan²⁸ with genome-wide variant associations derived from the same samples, binary outcomes and covariates. The resulting z-scores and effect sizes were compared with Pearson’s correlation to demonstrate the limitations of S-snTWAS in estimating accurate association effect sizes.

I-snTWAS and snGReX-PheWAS heterogeneity statistics were calculated for S-snTWAS significant (FDR-adjusted two-sided P ≤ 0.05) GTA effect sizes in R using the metafor package⁹⁵. Heterogeneity³⁰I²_{cell type} was calculated for genes imputable in at least two cell-type TIMs.

Targeted sn-eQTL analysis

To validate heterogeneity found in the same gene across cell types, we utilized PsychAD eQTLs calculated among EUR using multivariate multiple quantitative trait loci¹². We used the qtlPlots R package (v0.0.7) to visualize eQTLs⁹⁶. The top SNP for each cell type was extracted and linkage disequilibrium between them was examined using LDlink³¹.

We used joint effect on phenotype of eQTLs associated with a gene in mixed cohorts-pathways (JEPEGMIX2-P)³² to perform linkage disequilibrium-aware competitive pathway enrichment analysis. For this analysis, we used GRCh37-aligned GWAS summary statistics (prepared as above in ‘Enrichment of novel genes’) and snTIMs converted into JEPEGMIX2-P-compatible annotation files. In addition, we incorporated biological pathways from Gene Ontology (Gene Ontology 2015), accessed through Molecular Signatures Database (MSigDB) 5.1. JEPEGMIX2-P is designed to perform pathway enrichment analysis while accounting for linkage disequilibrium structures among genetic variants. For this study, the software was enhanced to conduct competitive analyses using the correlation adjusted mean rank (CAMERA)⁹⁷ gene set test procedure, allowing for a more refined understanding of gene–pathway associations; the publicly available binary executable file was updated to include this enhancement. The derived two-sided P values were FDR-adjusted²¹ among all pathways across all ancestry-specific TIMs, and an 0.05 FDR-adjusted two-sided P value threshold was used to determine significance. To obtain more conservative estimates of the number of significant pathways per cell type while maintaining specificity, we removed all significant pathways that were ‘parents’ of other significant pathways (referred to as pathway pruning in the main text). Finally, we note that we removed the MAPT locus (defined as chromosome 17: 44928498−56807609 in GRCh38) due to the presence of haplotypes that may bias results⁹⁸. Of note, these haplotypes have been found to have a large role in AD, PD and other disorders⁹⁸.

Comparison of MG AD snTWAS with differential gene expression signatures associated with AD pathology

We utilized transcriptional profiling data from freshly isolated primary human MG, along with the corresponding AD-relevant differentially expressed gene signatures and hierarchical co-expression networks¹⁹, which were constructed using multiscale embedded gene co-expression network analysis (MEGENA)⁹⁹. The MEGENA network was constructed from FACS-sorted CD45⁺ primary human MG from 189 autopsy samples (both with and without AD-associated pathology; comprising a partially intersecting set of the FACS-MG cohort and utilizing standard settings, including: Pearson’s correlation, un-signed and minimum module size of at least ten genes. For further analyses, we retained 306 co-expression modules with 50 or more genes each, and evaluated their enrichment for differentially expressed genes and S-snTWAS signatures using the competitive gene-set test CAMERA⁹⁷ implemented in the limma R package¹⁰⁰ with default settings, using pre-ranked z-scores from the differentially expressed gene and S-snTWAS analyses. The clinical dementia rating (CDR) phenotype indicates CDR ≥ 1 versus CDR = 0, and ‘Braak’ is a categorical phenotype separating high values ≥ 5 and low values ≤ 2. β-Amyloid density indicates the mean of this metric across five measured brain regions. Functional annotation enrichment of each MEGENA module using signatures from MSigDB (v7.2)¹⁰¹ was assessed using Fisher’s exact test. The relationship between differentially expressed genes and S-snTWAS was evaluated via Pearson’s correlations between the significance estimates (signed −log₁₀(two-sided P value)) of the enrichment analyses across all modules. The modules M406, M824 and M947 were the only modules with FDR-adjusted two-sided P ≤ 0.05 in the S-snTWAS enrichments via CAMERA pre-ranked (CameraPR).

To assess and visualize the concordance of z-scores between two datasets (for example, TWAS or summary statistics from competitive pathway enrichment analysis) at increasing significance thresholds, we performed the following analysis, which we termed PTCA. First, we scaled z-scores among all values in each dataset (using R’s base scale function, to normalize the values with a mean of 0 and a standard deviation of 1). Next, we matched identifiers (for example, trait–cell-type–gene or trait–cell-type–pathway) between the two datasets and restricted the final dataset to elements common to both datasets. Then, we performed a fixed-effect inverse-variance weighted meta-analysis using the z-scores of the dataset (standard error is assumed to be 1), and sorted the table based on the meta-analysis two-sided P value (increasing order). Finally, we measured the correlation between ordered and scaled z-scores in a step-wise restrictive manner (assuming step size = 10 and the dataset consists of 1,000 common elements, we assessed the correlation between the top 1,000 elements, followed by the top 990 elements, and so on).

Cross-disorder analysis

To compare snTWAS results across disorders, we used PTCA (described above) to assess how cross-trait z-score correlations strengthen as significance thresholds tighten, and two-sided Fisher’s exact tests to test enrichment of shared significant associations. For validation, we repeated the cross-disorder analyses with I-snTWAS in non-overlapping samples (no known overlap except for MDD) and required S-snTWAS-significant gene–cell-type associations to show concordant directions in I-snTWAS. This validation was not performed for trait pairs lacking adequate MVP power or for MDD, whose GWAS includes MVP participants. For all gene-level and pathway-level comparisons, we excluded the major histocompatibility complex and MAPT (GRCh38 chromosome 17: 44928498–56807609) regions due to haplotype complexity that can bias results (see ‘TWAS pathway enrichment analysis’)⁹⁸.

The data core team from the MVP handled genotyping, SNP imputation and initial quality control of genotypes. For this study, we utilized the MVP Release 4, which includes genotypes from 662,681 individuals. DNA was extracted from whole blood and genotyping was performed with the MVP 1.0 custom Axiom Array⁶⁰. The MVP data core called genotypes using assay for transposase-accessible chromatin (ATAC) ATAC Primer Tool (v2.11.3), Affy2vcf, PLINK⁹¹ and MVP 1.0 array library r6. Genotyping was validated using three plates of 1000 Genome samples. Rigorous genotype quality control such as plate normalization was used to improve the accuracy of genotype calling. Genotype imputation was performed using the TOPMed⁶¹ reference, SHAPEIT4 (v4.1.3)¹⁰² and Minimac4 (ref. ⁶²), and genetic principal components were generated using EIGENSOFT (v6)^103,104. We then performed sample-level and variant-level quality control on genotypes. Our quality control pipeline is based on previously established methods for other genetic analyses and optimized for maximizing input information for GReX estimation^6,74,105,106. We grouped samples by ancestry based on harmonized ancestry and race/ethnicity¹⁰⁷ (a method that integrates and harmonizes self-identified ancestry and genetic ancestry) into three ancestries (438,582 EUR, 112,346 AFR and 48,726 AMR). Using only autosomal variants, we filtered variants for MAF (0.01 or more), Hardy–Weinberg equilibrium P value (less than 1 × 10⁻⁶) and removed high disequilibrium regions (GRCh37; chromosome 6: 25500000−33500000, chromosome 8: 8135000−12000000 and chromosome 17: 40900000−45000000)¹⁰⁸. We then filtered for excess heterozygosity (four or more standard deviations from the mean), ambiguous sex (using both PLINK’s check-sex option, as well as sex chromosome aneuploidies assignments derived from plate intensity measurements) and relatedness (kinship-based inference for GWAS (KING) filter = 0.0884)^91,109. Finally, we limited variants to biallelic SNPs for individual imputation of GReX.

Phenotypes from International Classification of Diseases-9/10 codes were transformed into phecodes using Phecode Map (v1.2)^110,111. Individuals with at least two phecodes for a given trait were considered a case for both I-snTWAS and snGReX-PheWAS analyses. Owing to the enrichment of neuropsychiatric disorders within the MVP and their high comorbidity, individuals with fewer than two phecodes were considered a control for both analyses. For the I-snTWAS analysis, only traits with at least 2% case prevalence in EUR, AFR and AMR were utilized for power considerations (Supplementary Table 12), and AD was mapped to the phecode for ‘delirium dementia and amnestic disorders’, 290, due to the way AD is routinely coded in the Veteran Affairs medical system¹¹². For the snGReX-PheWAS analysis, traits with at least 500 cases and 500 controls were considered.

We estimated GReX in the MVP on an individual level using TIM-derived SNP predictor weights. Missing genotypes were replaced with double the in-sample MAF (representing the in-sample average genotype), rounded to the nearest integer.

Logistic regression analysis was performed while adjusting for sex, age and the top ten genetic principal components. Multiple test correction with FDR-adjusted two-sided P value²¹ (0.05 or less) was performed across all cell types, ancestries and traits. For comparison of I-snTWAS with S-snTWAS, we also utilized GWAS summary statistics for anxiety⁷⁶ and PTSD⁷⁵ (Fig. 4c and Supplementary Table 12).

Gene selection for snGReX-PheWAS

Owing to computational resource limitations, PheWAS analysis was performed on the top-ranked genes from the S-snTWAS and I-snTWAS. We selected significant genes from our S-snTWAS and I-snTWAS (FDR-adjusted two-sided P ≤ 0.05), and performed PheWAS on the top associations (Bonferroni-adjusted two-sided P ≤ 0.05 across all confidently imputable associations within each cell type in S-snTWAS or I-snTWAS and at least FDR-adjusted two-sided P ≤ 0.05 in the complementary snTWAS analysis; 623 unique genes) to maximize analytical yield while remaining within computational constraints. S-snTWAS genes were considered across all cell types in EUR. I-snTWAS genes were considered across all cell types and ancestries. This resulted in a total of 9,611 PheWASs (approximately 6.42% of all imputable coding gene–cell-type combinations) across 623 unique genes.

As in I-snTWAS, logistic regression analysis was performed while adjusting for sex, age and the top ten genetic principal components. FDR-adjustment was performed across all considered phecodes, 623 selected genes, 3 ancestries and all cell types for considered genes. Associations with FDR-adjusted two-sided P ≤ 0.05 (ref. ²¹) were considered significant.

To explore the pleiotropic effects on multiple phenotypes of the top snTWAS GTAs, we performed clustering of snGReX-PheWAS results. For interpretation of associations between brain disorders, we limited phenotype categories to ‘mental disorders’ and ‘neurological’. FDR adjustment was performed across all phecodes in these categories, as well as across all 623 PheWAS genes (in all cell types and ancestries) as typical. For every trait, we extracted the top 20 S-snTWAS and I-snTWAS significant genes (FDR-adjusted two-sided P ≤ 0.05 (ref. ²¹)) in association with the target phecode (Supplementary Table 12). We extracted only the top associated gene–cell-type combination for each gene to prevent clustering due to homogeneity of GReX across cell types. We then extracted up to the top 20 phecodes significantly associated (FDR-adjusted two-sided P ≤ 0.05 (ref. ²¹)) with these gene–cell-type combinations. Ward’s hierarchical agglomerative clustering¹¹³ was performed on the association z-scores with an implementation that preserves Ward’s criterion¹¹⁴ (ward.D2 method in R). Finally, we note that we removed the MAPT locus (defined as chromosome 17: 44928498−56807609 in GRCh38) due to the presence of haplotypes that may bias results (see ‘TWAS pathway enrichment analysis’)⁹⁸.

Abundance analysis of imputable genes

We hypothesized that genes uniquely identified by snBulk versus other snTIMs were lower in abundance. To investigate this, we subsetted imputable genes into those uniquely identified by snBulk (571) and those identified by at least one other snTIM (8,388). We calculated the mean number of transcripts for every gene among all EUR individuals with available expression data for the gene. To establish the statistical significance of the difference in mean transcript count, we used the Wilcoxon rank-sum test.

Out-of-sample validation of PsychAD snTIMs

Validation of class-level PsychAD snTIMs in ROSMAP

We performed external validation of EUR snTIMs using WGS and snRNA-seq from 396 ROSMAP¹⁷ participants of EUR. WGS aligned to GRCh37 was lifted over to GRCh38 (ref. ⁹³). To determine sample ancestry, the ROSMAP genotypes underwent standard genome-wide quality control procedures^115,116. We performed an initial round of variant quality control by removing variants with less than 0.95 call rate, to avoid biases in our sample quality control. Then, we proceeded with the removal of samples that fit any of the following criteria: call rate < 0.98, absolute value of the inbreeding coefficient < 0.2, genomic sex discrepancy with reported sex and formation of pairs with relatedness (PLINK’s calculated via identity by descent) > 0.1875 (ref. ⁹¹). We applied variant quality control, excluding markers with call rate < 0.98, and Hardy–Weinberg equilibrium P < 10⁻⁶. For ancestry validation and assignment, the resulting dataset was subsequently merged with the 1000 Genomes Project reference panel, keeping variants with MAF > 0.01 and pruned with a window size of 100, a step size of 5 and a pairwise R² threshold of 0.1. The samples were subsequently projected onto the eigenvectors computed on the 1000 Genomes Project samples using EIGENSOFT’s smartpca¹⁰³. Finally, we utilized quadratic discriminant analysis-determined EUR ROSMAP samples to impute individual-level GReX using PsychAD snTIMs (the ROSMAP WGS covered 99.7% of SNPs used by the PsychAD snTIMs).

Cells from the ROSMAP snRNA-seq dataset were annotated using the PsychAD taxonomy. We first subsampled 1 million cells from both the PsychAD and the ROSMAP snRNA-seq datasets, and then used the mapping between the two datasets to infer the labels for the full dataset. We used single-cell annotation using variational inference¹¹⁷ to perform reference-based label transfer. After subsetting the genes shared between two datasets, we used the scvi-tools package¹¹⁸ to train the scVI model based on the reference dataset and inferred cell-type labels (for example, class and subclass) on the query dataset. Models were run with 5 hidden layers and 10 latent variables, and the single-cell annotation using variational inference model was trained for 20 epochs with a minimal sample of 100 cells per cluster per epoch. Last, a transfer model was trained for 100 epochs and applied to query data to assign labels based on those the model was trained on from the reference. Label transfer achieved 91.6% accuracy when evaluated using true labels from the reference dataset. ROSMAP pseudobulk data were processed identically to PsychAD.

snRNA-seq profiles were clustered according to the PsychAD taxonomy and pseudobulked per cell type. Gene expression quality control matched the PsychAD pipeline above, and pseudobulk expression was PEER residualized. For each gene–cell-type combination present in both datasets, we computed the Pearson’s correlation across individuals between observed expression and GReX, defining out-of-sample R² as r². Agreement with internal cross-validated performance was assessed by Spearman’s correlation of out-of-sample R² versus R²_CV.

Validation of PsychAD subclass-MG in FACS-MG

To perform out-of-sample validation in FACS-MG, individual genotypes (filtering detailed in ‘S-snTWAS validation’ below) were leveraged to calculate GReX using the PsychAD EUR subclass-MG snTIM. GReX was then compared with PEER-residualized FACS-MG RNA-seq (see ‘S-snTWAS validation’ below). For each gene present in both datasets, we computed the Pearson’s correlation across individuals between observed expression and GReX, defining out-of-sample R² as r². Agreement with internal cross-validated performance was assessed by Spearman’s correlation of out-of-sample R² versus R²_CV.

We validated S-snTWAS associations using two external datasets: FACS-MG and Zeng-2024 (ref. ⁵¹). To prepare the FACS-MG TIM, we utilized similar protocols to the PsychAD snTIMs. FACS-MG cohort samples were derived from one of three different sub-cohort (based on the site of sample preparation). Genotypes were TOPMed-imputed and a EUR SNP reference panel similar to that used for PsychAD snTIM creation was used to initially select SNPs. The 1000 Genomes Project overlap was not implemented for the reference panel used in FACS-MG TIM creation due to the application of the TIM to a single GWAS: SNP utilization for the FACS-MG AD S-TWAS was more than 95%. Samples were filtered for missingness (0.01 or less) and relatedness (KING filter = 0.0884). To perform population stratification, variants were filtered for MAF (0.05 or more), missingness (0.01 or less) and Hardy–Weinberg equilibrium P value (less than 10⁻¹⁰). Variants were pruned using PLINK’s ‘--indep-pairwise’ function (1000 10 0.02). Twenty principal components were calculated and principal component analysis was used to determine EUR individuals (EUR selection ellipsoid defined using three standard deviations and three principal components). After ancestry filtering, sample-level quality control was performed. Variants were filtered for missingness (0.01 or less) before sample filtering for missingness (0.01 or less). Subsequently, only autosomal variants were retained, and variants were filtered for MAF (0.05 or more) and Hardy–Weinberg equilibrium P value (less than 10⁻⁶). High linkage disequilibrium regions were then filtered out (Supplementary Table 42) before pruning using PLINK’s ‘--indep-pairwise’ function (50 5 0.02). We then filtered for excess heterozygosity (three or more standard deviations from the mean) before filtering samples for relatedness (KING filter = 0.0884). After retaining only quality-controlled EUR samples, downstream variant filtering was performed identically to PsychAD. RNA-seq was processed identically to PsychAD snRNA-seq, including the two step PEER factor optimization on the three sub-cohorts. After filtering, we retained 271 EUR individuals with genotypes and RNA-seq. Finally, the FACS-MG TIM was filtered identically to PsychAD snTIMs.

FACS-MG TWAS was compared with EUR subclass-MG snTWAS using publicly available AD GWAS⁷⁸. For comparison parity, FDR adjustment was applied to all confidently imputable coding genes (after removing the major histocompatibility complex locus) in FACS-MG and subclass-MG separately. We then used z-scores from each to compare TWAS using Pearson’s correlation.

PsychAD class-level snTIMs were compared to Zeng-2024 snTIMs (Supplementary Table 39) using the same MDD GWAS summary statistics¹¹⁹. GWAS summary statistics were imputed for missing SNPs to ensure adequate coverage of PsychAD snTIM SNPs. Zeng-2024 TWAS was publicly available, and identical snTIM filtering to PsychAD was applied (R²_CV ≥ 0.01, FDR-adjusted two-sided P_CV ≤ 0.05). For comparison parity, FDR adjustment was applied to all confidently imputable coding genes (after removing the major histocompatibility complex locus) in Zeng-2024 and PsychAD separately. We then used z-scores from each to compare TWAS using Pearson’s correlation.

Genetic correlation analysis

We performed bivariate heritability analysis using linkage disequilibrium score regression^120,121 to assess the genetic correlation between the aforementioned seven NPD, four NDD and three SUD GWAS summary statistics. The datasets were munged to match the HapMap3 SNP allelic information, and the linkage disequilibrium weights were pre-calculated on the 1000 Genomes Project EUR dataset¹²².

Statistics and reproducibility

To demonstrate reproducibility, we performed multiple independent and internal validation analyses. Here we summarize five key examples. (1) We performed out-of-sample replication of PsychAD snTIMs in the ROSMAP and FACS-MG cohorts (Supplementary Fig. 26a,c). We correlated imputed GReX with observed PEER-residualized snRNA-seq expression and found strong concordance between out-of-sample performance (R²_out) and cross-validation performance (R²_CV; Spearman’s ρ = 0.63; N = 70,156; P = 3.91 × 10^−7,844, and ρ = 0.55; N = 1,853, P = 1.67 × 10⁻¹⁴⁹, respectively). (2) We compared S-snTWAS results derived from PsychAD snTIMs with those obtained using independent transcriptomic imputation models, including previously published snTIMs⁵¹ and a MG-specific TIM derived from FACS-MG¹⁹ (Supplementary Fig. 26b,d). In both cases, we observed strong concordance (Pearson’s r_(13,760) = 0.78; P = 1.34 × 10^−2,772 and Pearson’s r₍₇₃₃₎ = 0.87; P = 1.32 × 10⁻²³¹, respectively). (3) To address potential bias in the cross-disorder analysis from sample overlap across GWAS datasets, we replicated the analysis replacing S-snTWAS data with independent MVP I-snTWAS data and evaluated sign concordance of results (Supplementary Fig. 14). Overall, 90% of cross-disorder associations (5,162 of 5,721) showed consistent direction of effect. (4) We further compared S-snTWAS with MVP EUR I-snTWAS (Fig. 5a). Although the overall Pearson correlation across all associations was modest, PTCA demonstrated strong concordance among top-ranked associations. (5) Finally, snTIM training incorporated fivefold cross-validation, consistent with the broader use of cross-validation for evaluating PrediXcan-style transcriptomic imputation models.

Reporting summary

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

Data availability

All results are included either in the main text or provided in Supplementary Tables or external data²⁰. We are also making snTIMs and accompanying SNP covariance matrices publicly available on Zenodo¹²³ (https://doi.org/10.5281/zenodo.17435917) and Synapse²⁰ (https://doi.org/10.7303/syn63181047). The PsychAD single-nucleus transcriptomics data are available via the AD Knowledge Portal (https://doi.org/10.7303/9618136). Source data are provided with this paper.

Code availability

This project utilized publicly available code as described in the Reporting Summary.

Brainstorm Consortium et al. Analysis of shared heritability in common disorders of the brain. Science 360, eaap8757 (2018).
Girdhar, K. et al. Cell-specific histone modification maps in the human frontal lobe link schizophrenia risk to the neuronal epigenome. Nat. Neurosci. 21, 1126–1136 (2018).
Girdhar, K. et al. Chromatin domain alterations linked to 3D genome organization in a large cohort of schizophrenia and bipolar disorder brains. Nat. Neurosci. 25, 474–483 (2022).
Hauberg, M. E. et al. Common schizophrenia risk variants are enriched in open chromatin regions of human glutamatergic neurons. Nat. Commun. 11, 5581 (2020).
Kosoy, R. et al. Genetics of the human microglia regulome refines Alzheimer’s disease risk loci. Nat. Genet. 54, 1145–1154 (2022).
Voloudakis, G. et al. A translational genomics approach identifies IL10RB as the top candidate gene target for COVID-19 susceptibility. NPJ Genom. Med. 7, 52 (2022).
Gamazon, E. R. et al. A gene-based association method for mapping traits using reference transcriptome data. Nat. Genet. 47, 1091–1098 (2015).
Gusev, A. et al. Integrative approaches for large-scale transcriptome-wide association studies. Nat. Genet. 48, 245–252 (2016).
Zhang, W. et al. Integrative transcriptome imputation reveals tissue-specific and shared biological mechanisms mediating susceptibility to complex traits. Nat. Commun. 10, 3834 (2019).
Wainberg, M. et al. Opportunities and challenges for transcriptome-wide association studies. Nat. Genet. 51, 592–599 (2019).
Lee, D. et al. Single-cell atlas of transcriptomic vulnerability across brain disorders. Nature https://doi.org/10.1038/s41586-025-09573-z (2026).
Zeng, B. et al. Single-nucleus atlas of cell-type specific genetic regulation in the human brain. Nat. Genet. https://doi.org/10.1038/s41588-026-02733-5 (2026).
Zeng, B. et al. Genetic regulation of cell type-specific chromatin accessibility shapes brain disease etiology. Science 384, eadh4265 (2024).
Fullard, J. F. et al. An atlas of chromatin accessibility in the adult human brain. Genome Res. 28, 1243–1252 (2018).
Fullard, J. F. et al. Population-scale cross-disorder atlas of the human prefrontal cortex at single-cell resolution. Sci. Data 12, 954 (2025).
Fullard, J. F. et al. Single-nucleus transcriptome analysis of human brain immune response in patients with severe COVID-19. Genome Med. 13, 118 (2021).
Mathys, H. et al. Single-cell atlas reveals correlates of high cognitive function, dementia, and resilience to Alzheimer’s disease pathology. Cell 186, 4365–4385.e27 (2023).
Humphrey, J. et al. Long-read RNA sequencing atlas of human microglia isoforms elucidates disease-associated genetic regulation of splicing. Nat. Genet. 57, 604–615 (2025).
Kosoy, R. et al. Alzheimer’s disease transcriptional landscape in ex vivo human microglia. Nat. Neurosci. 28, 1830–1843 (2025).
Venkatesh, S. PsychAD_snTWAS. Preprint at Synapse https://doi.org/10.7303/syn63181047 (2024).
Benjamini, Y., Drai, D., Elmer, G., Kafkafi, N. & Golani, I. Controlling the false discovery rate in behavior genetics research. Behav. Brain Res. 125, 279–284 (2001).
McKusick, V. A. Mendelian Inheritance in Man: a Catalog of Human Genes and Genetic Disorders (JHU Press, 1998).
Watanabe, K., Taskesen, E., van Bochoven, A. & Posthuma, D. Functional mapping and annotation of genetic associations with FUMA. Nat. Commun. 8, 1826 (2017).
de Leeuw, C. A., Mooij, J. M., Heskes, T. & Posthuma, D. MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput. Biol. 11, e1004219 (2015).
Mancuso, N. et al. Probabilistic fine-mapping of transcriptome-wide association studies. Nat. Genet. 51, 675–682 (2019).
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).
Baker, M. R., Lee, A. S. & Rajadhyaksha, A. M. L-type calcium channels and neuropsychiatric diseases: insights into genetic risk variant-associated genomic regulation and impact on brain development. Channels 17, 2176984 (2023).
Barbeira, A. N. et al. Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. Nat. Commun. 9, 1825 (2018).
Yoon, S., Baik, B., Park, T. & Nam, D. Powerful p-value combination methods to detect incomplete association. Sci. Rep. 11, 6980 (2021).
Higgins, J. P. T. & Thompson, S. G. Quantifying heterogeneity in a meta-analysis. Stat. Med. 21, 1539–1558 (2002).
Machiela, M. J. & Chanock, S. J. LDlink: a web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants. Bioinformatics 31, 3555–3557 (2015).
Chatzinakos, C. et al. TWAS pathway method greatly enhances the number of leads for uncovering the molecular underpinnings of psychiatric disorders. Am. J. Med. Genet. B Neuropsychiatr. Genet. 183, 454–463 (2020).
Jimerson, D. C. et al. CSF calcium: clinical correlates in affective illness and schizophrenia. Biol. Psychiatry 14, 37–51 (1979).
Lanni, C. et al. Zyxin is a novel target for β-amyloid peptide: characterization of its role in Alzheimer’s pathogenesis. J. Neurochem. 125, 790–799 (2013).
Xie, C. et al. A shared neural basis underlying psychiatric comorbidity. Nat. Med. 29, 1232–1242 (2023).
Hesdorffer, D. C. Comorbidity between neurological illness and psychiatric disorders. CNS Spectr. 21, 230–238 (2016).
Reynolds, R. H. et al. Local genetic correlations exist among neurodegenerative and neuropsychiatric diseases. NPJ Parkinsons Dis. 9, 70 (2023).
Bhattacharya, A. et al. Best practices for multi-ancestry, meta-analytic transcriptome-wide association studies: lessons from the Global Biobank Meta-analysis Initiative. Cell Genom. 2, 100180 (2022).
Keys, K. L. et al. On the cross-population generalizability of gene expression prediction models. PLoS Genet. 16, e1008927 (2020).
Clarke, S. L. et al. Race and ethnicity stratification for polygenic risk score analyses may mask disparities in Hispanics. Circulation 146, 265–267 (2022).
Langhammer, F. et al. Deregulated ion channels contribute to RHOBTB2-associated developmental and epileptic encephalopathy. Hum. Mol. Genet. 34, 639–650 (2025).
O’Connell, K. S. et al. Genomics yields biological and phenotypic insights into bipolar disorder. Nature 639, 968–975 (2025).
Liang, Y. et al. Polygenic transcriptome risk scores (PTRS) can improve portability of polygenic risk scores across ancestries. Genome Biol. 23, 23 (2022).
Hasin, D. S., Stinson, F. S., Ogburn, E. & Grant, B. F. Prevalence, correlates, disability, and comorbidity of DSM-IV alcohol abuse and dependence in the United States: results from the National Epidemiologic Survey on Alcohol and Related Conditions. Arch. Gen. Psychiatry 64, 830–842 (2007).
Hunt, G. E., Malhi, G. S., Cleary, M., Lai, H. M. X. & Sitharthan, T. Comorbidity of bipolar and substance use disorders in national surveys of general populations, 1990-2015: systematic review and meta-analysis. J. Affect. Disord. 206, 321–330 (2016).
Schwarzinger, M. et al. Contribution of alcohol use disorders to the burden of dementia in France 2008-13: a nationwide retrospective cohort study. Lancet Public Health 3, e124–e132 (2018).
de Paiva Lopes, K. et al. Genetic analysis of the human microglial transcriptome across brain regions, aging and disease pathologies. Nat. Genet. 54, 4–17 (2022).
Young, A. M. H. et al. A map of transcriptional heterogeneity and regulatory variation in human microglia. Nat. Genet. 53, 861–868 (2021).
Jaffe, A. E. et al. Profiling gene expression in the human dentate gyrus granule cell layer reveals insights into schizophrenia and its genetic risk. Nat. Neurosci. 23, 510–519 (2020).
Bryois, J. et al. Cell-type-specific cis-eQTLs in eight human brain cell types identify novel risk genes for psychiatric and neurological disorders. Nat. Neurosci. 25, 1104–1112 (2022).
Zeng, L. et al. A single-nucleus transcriptome-wide association study implicates novel genes in depression pathogenesis. Biol. Psychiatry 96, 34–43 (2024).
Patro, N. & Patro, I. in The Biology of Glial Cells: Recent Advances (eds Patro, I. et al.) 143–170 (Springer Singapore, 2022).
Andersen, M. S. et al. Heritability enrichment implicates microglia in Parkinson’s disease pathogenesis. Ann. Neurol. 89, 942–951 (2021).
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).
Mogil, L. S. et al. Genetic architecture of gene expression traits across diverse populations. PLoS Genet. 14, e1007586 (2018).
Ping, J. et al. Using genome and transcriptome data from African-ancestry female participants to identify putative breast cancer susceptibility genes. Nat. Commun. 15, 3718 (2024).
Wen, C. et al. Cross-ancestry atlas of gene, isoform, and splicing regulation in the developing human brain. Science 384, eadh0829 (2024).
1000 Genomes Project Consortium et al. A global reference for human genetic variation. Nature 526, 68–74 (2015).
Fairley, S., Lowy-Gallego, E., Perry, E. & Flicek, P. The International Genome Sample Resource (IGSR) collection of open human genomic variation resources. Nucleic Acids Res. 48, D941–D947 (2020).
Gaziano, J. M. et al. Million Veteran Program: a mega-biobank to study genetic influences on health and disease. J. Clin. Epidemiol. 70, 214–223 (2016).
Taliun, D. et al. Sequencing of 53,831 diverse genomes from the NHLBI TOPMed Program. Nature 590, 290–299 (2021).
Das, S. et al. Next-generation genotype imputation service and methods. Nat. Genet. 48, 1284–1287 (2016).
Fuchsberger, C., Abecasis, G. R. & Hinds, D. A. minimac2: Faster genotype imputation. Bioinformatics 31, 782–784 (2015).
Phan, L. et al. ALFA: allele frequency aggregator. NCBI https://www.ncbi.nlm.nih.gov/snp/docs/gsr/alfa/#citing-this-project (2020).
Hoffman, G. E. et al. Efficient differential expression analysis of large-scale single-cell transcriptomics data using Dreamlet. Nat. Commun. https://doi.org/10.1038/s41467-026-75680-8 (2026).
Cunningham, F. et al. Ensembl 2022. Nucleic Acids Res. 50, D988–D995 (2022).
Nassar, L. R. et al. The UCSC Genome Browser database: 2023 update. Nucleic Acids Res. 51, D1188–D1195 (2023).
Stegle, O., Parts, L., Piipari, M., Winn, J. & Durbin, R. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses. Nat. Protoc. 7, 500–507 (2012).
Ongen, H., Buil, A., Brown, A. A., Dermitzakis, E. T. & Delaneau, O. Fast and efficient QTL mapper for thousands of molecular phenotypes. Bioinformatics 32, 1479–1485 (2016).
Nyasimi, F. hakyimlab/PredictDB-Tutorial: tutorial for running the PredictDB pipeline with GEUVADIS data. GitHub https://github.com/hakyimlab/PredictDB-Tutorial (2021).
Mancuso, N. et al. Large-scale transcriptome-wide association study identifies new prostate cancer risk regions. Nat. Commun. 9, 4079 (2018).
Wilkinson, G. N. & Rogers, C. E. Symbolic description of factorial models for analysis of variance. J. R. Stat. Soc. C Appl. Stat. 22, 392–399 (1973).
Chambers, J. M. & Hastie, T. Statistical Models in S (Chapman & Hall/CRC, 1992).
Burstein, D. et al. Genome-wide analysis of a model-derived binge eating disorder phenotype identifies risk loci and implicates iron metabolism. Nat. Genet. 55, 1462–1470 (2023).
Nievergelt, C. M. et al. International meta-analysis of PTSD genome-wide association studies identifies sex- and ancestry-specific genetic risk loci. Nat. Commun. 10, 4558 (2019).
Dönertaş, H. M., Fabian, D. K., Valenzuela, M. F., Partridge, L. & Thornton, J. M. Common genetic associations between age-related diseases. Nat. Aging 1, 400–412 (2021).
International Obsessive Compulsive Disorder Foundation Genetics Collaborative (IOCDF-GC) and OCD Collaborative Genetics Association Studies (OCGAS). Revealing the complex genetic architecture of obsessive-compulsive disorder using meta-analysis. Mol. Psychiatry 23, 1181–1188 (2018).
Bellenguez, C. et al. New insights into the genetic etiology of Alzheimer’s disease and related dementias. Nat. Genet. 54, 412–436 (2022).
van Rheenen, W. et al. Common and rare variant association analyses in amyotrophic lateral sclerosis identify 15 risk loci with distinct genetic architectures and neuron-specific biology. Nat. Genet. 53, 1636–1648 (2021).
Watson, H. J. et al. Genome-wide association study identifies eight risk loci and implicates metabo-psychiatric origins for anorexia nervosa. Nat. Genet. 51, 1207–1214 (2019).
Demontis, D. et al. Genome-wide analyses of ADHD identify 27 risk loci, refine the genetic architecture and implicate several cognitive domains. Nat. Genet. 55, 198–208 (2023).
Mullins, N. et al. Genome-wide association study of more than 40,000 bipolar disorder cases provides new insights into the underlying biology. Nat. Genet. 53, 817–829 (2021).
Jansen, P. R. et al. Genome-wide analysis of insomnia in 1,331,010 individuals identifies new risk loci and functional pathways. Nat. Genet. 51, 394–403 (2019).
Als, T. D. et al. Depression pathophysiology, risk prediction of recurrence and comorbid psychiatric disorders using genome-wide analyses. Nat. Med. 29, 1832–1844 (2023).
International Multiple Sclerosis Genetics Consortium. Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science 365, eaav7188 (2019).
Nalls, M. A. et al. Identification of novel risk loci, causal insights, and heritable risk for Parkinson’s disease: a meta-analysis of genome-wide association studies. Lancet Neurol. 18, 1091–1102 (2019).
Trubetskoy, V. et al. Mapping genomic loci implicates genes and synaptic biology in schizophrenia. Nature 604, 502–508 (2022).
Sanchez-Roige, S. et al. Genome-wide association study meta-analysis of the alcohol use disorders identification test (AUDIT) in two population-based cohorts. Am. J. Psychiatry 176, 107–118 (2019).
Murphy, A. E., Schilder, B. M. & Skene, N. G. MungeSumstats: a Bioconductor package for the standardization and quality control of many GWAS summary statistics. Bioinformatics 37, 4593–4596 (2021).
Pasaniuc, B. et al. Fast and accurate imputation of summary statistics enhances evidence of functional enrichment. Bioinformatics 30, 2906–2914 (2014).
Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).
Hoffman, G. GabrielHoffman/imputez: impute z-statistics for missing tests using observed z-statistics and the correlation matrix between z-statistics. GitHub https://github.com/GabrielHoffman/imputez (2023).
liftOver. Bioconductor http://bioconductor.org/packages/liftOver/ (2026).
Lu, Z. et al. Multi-ancestry fine-mapping improves precision to identify causal genes in transcriptome-wide association studies. Am. J. Hum. Genet. 109, 1388–1404 (2022).
Viechtbauer, W. Conducting meta-analyses in R with themetafor package. J. Stat. Softw. https://doi.org/10.18637/jss.v036.i03 (2010).
Hoffman, G. GabrielHoffman/qtlPlots: create QTL plots. GitHub https://github.com/GabrielHoffman/qtlPlots (2025).
Wu, D. & Smyth, G. K. Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res. 40, e133 (2012).
Pedicone, C., Weitzman, S. A., Renton, A. E. & Goate, A. M. Unraveling the complex role of MAPT-containing H1 and H2 haplotypes in neurodegenerative diseases. Mol. Neurodegener. 19, 43 (2024).
Song, W.-M. & Zhang, B. Multiscale embedded gene co-expression network analysis. PLoS Comput. Biol. 11, e1004574 (2015).
Smyth, G. K. in Bioinformatics and Computational Biology Solutions Using R and Bioconductor (eds. Gentleman, R. et al.) 397–420 (Springer, 2005).
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005).
Delaneau, O., Zagury, J.-F., Robinson, M. R., Marchini, J. L. & Dermitzakis, E. T. Accurate, scalable and integrative haplotype estimation. Nat. Commun. 10, 5436 (2019).
Patterson, N., Price, A. L. & Reich, D. Population structure and eigenanalysis. PLoS Genet. 2, e190 (2006).
Price, A. L. et al. Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909 (2006).
Marees, A. T. et al. A tutorial on conducting genome-wide association studies: quality control and statistical analysis. Int. J. Methods Psychiatr. Res. 27, e1608 (2018).
Voloudakis, G. et al. Neuropsychiatric polygenic scores are weak predictors of professional categories. Nat. Hum. Behav. https://doi.org/10.1038/s41562-024-02074-5 (2024).
Fang, H. et al. Harmonizing genetic ancestry and self-identified race/ethnicity in genome-wide association studies. Am. J. Hum. Genet. 105, 763–772 (2019).
Jiang, L. et al. Genome-wide association analyses of common infections in a large practice-based biobank. BMC Genomics 23, 672 (2022).
Manichaikul, A. et al. Robust relationship inference in genome-wide association studies. Bioinformatics 26, 2867–2873 (2010).
Wei, W.-Q. et al. Evaluating phecodes, clinical classification software, and ICD-9-CM codes for phenome-wide association studies in the electronic health record. PLoS ONE 12, e0175508 (2017).
Wu, P. et al. Mapping ICD-10 and ICD-10-CM codes to phecodes: workflow development and initial evaluation. JMIR Med. Inform. 7, e14325 (2019).
Sherva, R. et al. African ancestry GWAS of dementia in a large military cohort identifies significant risk loci. Mol. Psychiatry 28, 1293–1302 (2023).
Ward, J. H. Jr. Hierarchical grouping to optimize an objective function. J. Am. Stat. Assoc. 58, 236–244 (1963).
Murtagh, F. & Legendre, P. Ward’s hierarchical agglomerative clustering method: which algorithms implement Ward’s criterion? J. Classification 31, 274–295 (2014).
Truong, V. Q. et al. Quality control procedures for genome-wide association studies. Curr. Protoc. 2, e603 (2022).
Turner, S. et al. Quality control procedures for genome-wide association studies. Curr. Protoc. Hum. Genet. https://doi.org/10.1002/0471142905.hg0119s68 (2011).
Xu, C. et al. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models. Mol. Syst. Biol. 17, e9620 (2021).
Lopez, R., Regier, J., Cole, M. B., Jordan, M. I. & Yosef, N. Deep generative modeling for single-cell transcriptomics. Nat. Methods 15, 1053–1058 (2018).
Wray, N. R. et al. Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression. Nat. Genet. 50, 668–681 (2018).
Bulik-Sullivan, B. K. et al. LD score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat. Genet. 47, 291–295 (2015).
Bulik-Sullivan, B. et al. An atlas of genetic correlations across human diseases and traits. Nat. Genet. 47, 1236–1241 (2015).
Gazal, S. S-LDSC reference files. Zenodo https://doi.org/10.5281/zenodo.10515792 (2024).
Venkatesh, S. PsychAD single-nucleus transcriptomic imputation models [data set]. Zenodo https://doi.org/10.5281/zenodo.17435917 (2025).

Download references

We thank the MVP staff, researchers and volunteers, who have contributed to MVP, and especially those who previously served their country in the military and now generously agreed to enroll in the study (see https://mvp.va.gov for more information). The underlying work was based on data from the MVP, Office of Research and Development, Veterans Health Administration, and was supported by the Veterans Administration MVP award #000. The citation for MVP is Gaziano, J. M. et al.⁶⁰, and P.R. supervised this work. Study data from the ROSMAP study (https://doi.org/10.7303/9618239) were provided by the RADC, Rush University Medical Center, Chicago; additional phenotypic data can be requested (www.radc.rush.edu).

This study was supported by the US National Institutes of Health (NIH) under award numbers R01AG067025 (to P.R. and V.H.), R01AG082185 (to P.R., D.L. and V.H.), R01AG078657 (to G.V.), K08MH122911 (to G.V.), R01AG065582 (to P.R. and V.H.), R01MH125246 (to P.R.), R01AG050986 (to P.R.), R01AG095776 (to P.R.), U24AG087563 (to P.R.) and T32MH087004 (to K.T.). This work was supported in part through the computational and data resources and staff expertise provided by Scientific Computing and Data at the Icahn School of Medicine at Mount Sinai and supported by the Clinical and Translational Science Award grant UL1TR004419 from the National Center for Advancing Translational Sciences. Study data from The Mount Sinai PsychAD Study (https://doi.org/10.7303/9618136) were generated from the postmortem brain tissue provided by the Mount Sinai Brain Bank, the RADC (funding: P30AG10161, P30AG72975, R01AG15819, R01AG17917, R01AG22018, U01AG46152 and U01AG61356) and the NIMH-IRP Human Brain Collection Core (project #ZIC MH002903). The Mount Sinai Brain Bank specimens were provided through the NIH NeuroBioBank and supported by NIMH-75N95019C00049. Data collection was supported by NIA grant R01AG067025 and was conducted by the principal investigators P.R., V.H., S.F. and D.W. ROSMAP data collection was supported through funding by NIA grants P30AG10161 (ROS), R01AG15819 (ROSMAP; genomics and RNAseq), R01AG17917 (MAP), R01AG30146, R01AG36042 (5hC methylation and ATAC-seq), RC2AG036547 (H3K9Ac), R01AG36836 (RNA-seq), R01AG48015 (monocyte RNA-seq) RF1AG57473 (snRNA-seq), U01AG32984 (genomic and whole-exome sequencing), U01AG46152 (ROSMAP AMP-AD and targeted proteomics), U01AG46161 (TMT proteomics), U01AG61356 (WGS, targeted proteomics and ROSMAP AMP-AD), P30AG072975, the Illinois Department of Public Health (ROSMAP) and the Translational Genomics Research Institute (genomic). snRNA-seq data generation was funded by NIH grants U01AG061356, RF1AG057473 and U01AG046152 as part of the AMP-AD Consortium, as well as NIH grants R01AG066831 and U01AG072572.

Author information

Author notes

These authors contributed equally: Roman Kosoy, Zhenyi Wu
These authors jointly supervised this work: Georgios Voloudakis, Panos Roussos

Authors and Affiliations

Department of Psychiatry, Icahn School of Medicine at Mount Sinai, New York, NY, USA
Sanan Venkatesh, Roman Kosoy, Zhenyi Wu, Marios Anyfantakis, Christian Dillard, Prashant N. M., David Burstein, Deepika Mathur, Bukola Ajanaku, Fotis Tsetsos, Biao Zeng, Sonali Gupta, Rachel Bercovitch, Aram Hong, Clara Casey, Marcela Alvia, Zhiping Shao, Stathis Argyriou, Karen Therrien, Christian Porras, Collin Spencer, Hui Yang, Jaroslav Bendl, Jennifer Monteiro Fortes, Lyra Sheu, Marios Anyfantakis, Maxim Signaevsky, Mikaela Koutrouli, Milos Pjanic, Nicolas Y. Masse, Pavel Katsel, Pengfei Dong, Sanan Venkatesh, Sarah R. Murphy, Seon Kinrot, Steven P. Kleopoulos, Tereza Clarence, Vahram Haroutunian, Xinyi Wang, Zhenyi Wu, Vahram Haroutunian, Kiran Girdhar, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Georgios Voloudakis & Panos Roussos
Center for Disease Neurogenomics, Icahn School of Medicine at Mount Sinai, New York, NY, USA
Sanan Venkatesh, Roman Kosoy, Zhenyi Wu, Marios Anyfantakis, Christian Dillard, Prashant N. M., David Burstein, Deepika Mathur, Bukola Ajanaku, Fotis Tsetsos, Biao Zeng, Sonali Gupta, Rachel Bercovitch, Aram Hong, Clara Casey, Marcela Alvia, Zhiping Shao, Stathis Argyriou, Karen Therrien, Christian Porras, Collin Spencer, Hui Yang, Jaroslav Bendl, Jennifer Monteiro Fortes, Lyra Sheu, Marios Anyfantakis, Mikaela Koutrouli, Milos Pjanic, Nicolas Y. Masse, Pengfei Dong, Sanan Venkatesh, Sarah R. Murphy, Seon Kinrot, Steven P. Kleopoulos, Tereza Clarence, Xinyi Wang, Zhenyi Wu, Vahram Haroutunian, Kiran Girdhar, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Georgios Voloudakis & Panos Roussos
Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, NY, USA
Sanan Venkatesh, Roman Kosoy, Zhenyi Wu, Marios Anyfantakis, Christian Dillard, Prashant N. M., David Burstein, Deepika Mathur, Bukola Ajanaku, Fotis Tsetsos, Biao Zeng, Sonali Gupta, Rachel Bercovitch, Aram Hong, Clara Casey, Marcela Alvia, Zhiping Shao, Stathis Argyriou, Karen Therrien, Christian Porras, Collin Spencer, Fotios Tsetsos, Hui Yang, Jaroslav Bendl, Jennifer Monteiro Fortes, Lyra Sheu, Marios Anyfantakis, Maxim Signaevsky, Mikaela Koutrouli, Milos Pjanic, Nicolas Y. Masse, Pengfei Dong, Sanan Venkatesh, Sarah R. Murphy, Seon Kinrot, Steven P. Kleopoulos, Tereza Clarence, Vahram Haroutunian, Xinyi Wang, Zhenyi Wu, Vahram Haroutunian, Kiran Girdhar, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Georgios Voloudakis & Panos Roussos
Sanan Venkatesh, Roman Kosoy, Zhenyi Wu, Marios Anyfantakis, Christian Dillard, Prashant N. M., David Burstein, Deepika Mathur, Bukola Ajanaku, Fotis Tsetsos, Biao Zeng, Sonali Gupta, Rachel Bercovitch, Aram Hong, Clara Casey, Marcela Alvia, Zhiping Shao, Stathis Argyriou, Karen Therrien, Christian Porras, Collin Spencer, Fotios Tsetsos, Hui Yang, Jaroslav Bendl, Jennifer Monteiro Fortes, Lyra Sheu, Marios Anyfantakis, Mikaela Koutrouli, Milos Pjanic, Nicolas Y. Masse, Pengfei Dong, Sanan Venkatesh, Sarah R. Murphy, Seon Kinrot, Steven P. Kleopoulos, Tereza Clarence, Xinyi Wang, Zhenyi Wu, Vahram Haroutunian, Kiran Girdhar, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Georgios Voloudakis & Panos Roussos
Mental Illness Research, Education, and Clinical Center (VISN 2 South), James J. Peters VA Medical Center, Bronx, NY, USA
Sanan Venkatesh, Zhenyi Wu, Marios Anyfantakis, David Burstein, Deepika Mathur, Vahram Haroutunian, Vahram Haroutunian, Jaroslav Bendl, Gabriel E. Hoffman, Georgios Voloudakis & Panos Roussos
David Burstein, Deepika Mathur, Fotis Tsetsos, Fotios Tsetsos, Sanan Venkatesh, Gabriel E. Hoffman, Georgios Voloudakis & Panos Roussos
Institute for Genomics in Health (IGH), SUNY Downstate Health Sciences University, Brooklyn, NY, USA
Chris Chatzinakos & Tim Bigdeli
VA New York Harbor Healthcare System, Brooklyn, NY, USA
Department of Epidemiology and Biostatistics, School of Public Health, SUNY Downstate Health Sciences University, Brooklyn, NY, USA
David A. Bennett & Tim Bigdeli
Human Brain Collection Core, National Institute of Mental Health-Intramural Research Program, Bethesda, MD, USA
Pavan K. Auluck, Pavan Auluck & Stefano Marenco
Rush Alzheimer’s Disease Center, Rush University Medical Center, Chicago, IL, USA
David A. Bennett, Lisa L. Barnes & David A. Bennett
Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA
Athan Z. Li, Daifeng Wang, Gennadi Ryan, Noah Cohen Kalafut & Sayali A. Alatkar
Waisman Center, University of Wisconsin-Madison, Madison, WI, USA
Athan Z. Li, Chenfeng He, Chirag Gupta, Daifeng Wang, Jerome J. Choi, Kalpana H. Arachchilage, Noah Cohen Kalafut, Pramod B. Chandrashekar, Saniya Khullar, Sayali A. Alatkar, Ting Jin & Xiang Huang
Chenfeng He, Daifeng Wang, Kalpana H. Arachchilage, Pramod B. Chandrashekar, Saniya Khullar & Ting Jin
Novo Nordisk Foundation Center for Protein Research, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark
Chirag Gupta, Lars J. Jensen & Mikaela Koutrouli
Department of Psychiatry, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA
Colleen A. McClung & Madeline R. Scott
Taube/Koret Center for Neurodegenerative Disease Research, Gladstone Institutes, San Francisco, CA, USA
Gennadi Ryan, Monika Ahirwar & Steven Finkbeiner
Department of Population Health Sciences, University of Wisconsin-Madison, Madison, WI, USA
Jerome J. Choi
Department of Neurological Sciences, Rush University Medical Center, Chicago, IL, USA
Lisa L. Barnes
Vanderbilt Genetics Institute, Vanderbilt University Medical Center, Nashville, TN, USA
Logan C. Dumitrescu & Timothy J. Hohman
Vanderbilt Memory & Alzheimer’s Center, Vanderbilt University Medical Center, Nashville, TN, USA
Logan C. Dumitrescu & Timothy J. Hohman
Center for Systems and Therapeutics, Gladstone Institutes, San Francisco, CA, USA
Monika Ahirwar, Steven Finkbeiner & Vivek G. Ramaswamy
Department of Neurology, University of California San Francisco, San Francisco, CA, USA
Steven Finkbeiner & Vivek G. Ramaswamy
Department of Physiology, University of California San Francisco, San Francisco, CA, USA
Neuroscience and Biomedical Sciences Graduate Programs, University of California San Francisco, San Francisco, CA, USA
Department of Neuroscience, Icahn School of Medicine at Mount Sinai, New York, NY, USA

View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
Prashant N. M.
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
David A. Bennett
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:
John F. Fullard
View author publications
Search author on:
Gabriel E. Hoffman
View author publications
Search author on:
View author publications
Search author on:
View author publications
Search author on:

, Athan Z. Li
, Colleen A. McClung
, David A. Bennett
, Gabriel E. Hoffman
, Jerome J. Choi
, John F. Fullard
, Kalpana H. Arachchilage
, Lars J. Jensen
, Lisa L. Barnes
, Logan C. Dumitrescu
, Madeline R. Scott
, Nicolas Y. Masse
, Pavan K. Auluck
, Pramod B. Chandrashekar
, Prashant N. M.
, Sarah R. Murphy
, Sayali A. Alatkar
, Steven P. Kleopoulos
, Timothy J. Hohman
, Vivek G. Ramaswamy

S.V., G.V. and P.R. conceptualized the study and came up with the study design. S.V., R.K., C.D., P.N.M., D.B., D.M., C. Chatzinakos, B.A., B.Z., S.G., R.B., A.H., C. Casey, M. Anyfantakis, Z.S., S.A., K.T., T.B., P.A., D.A.B., S.M., V.H., K.G., J.B., D.L., J.F.F., G.E.H., G.V. and P.R. contributed data or analysis tools. S.V., R.K., Z.W., M. Alvia, F.T. and G.V. performed the analyses. S.V., J.F.F., G.V. and P.R. wrote the manuscript with input from all authors.

Corresponding authors

Correspondence to Georgios Voloudakis or Panos Roussos.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

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

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Overview of study design.

Schematic of the overall analysis workflow, from postmortem dorsolateral prefrontal cortex sampling and single-nucleus RNA sequencing (snRNA-seq) in PsychAD, through construction of ancestry- and cell-type-specific single-nucleus transcriptomic imputation models (snTIMs), to summary-level snTWAS across 12 brain disorder genome-wide association studies (GWAS) and individual-level single-nucleus genetically regulated gene expression phenome-wide association studies (snGReX-PheWAS) in the Million Veteran Program (MVP).

Extended Data Fig. 2 PsychAD cell-type hierarchy.

Schematic of the dorsolateral prefrontal cortex (DLPFC) cell-type hierarchy. Cell-types are organized into homogenate (Bulk/snBulk), Class, and Subclass levels. Single-nucleus transcriptomic imputation models (snTIMs) based on single-nucleus RNA-seq (snRNA-seq) of the PsychAD cohort are split into 3 levels depending on level of resolution (snBulk, Class and Subclass) whereas a previously trained TIM for the same region using homogenate RNA-seq is limited to the homogenate level (Bulk). Class-OPC, Astro, Endo, and Oligo are not further subdivided at the subclass level and are therefore included in both class and subclass aggregate analyses. This results in 27 non-overlapping cellular populations across 32 cellular populations in the PsychAD cohort. Full names of all cell-types are defined in Supplementary Table 3.

Extended Data Fig. 3 Cell-type-specific genes and pathways in single-nucleus transcriptome-wide association studies (snTWAS).

A, Heatmap of multivariate adaptive shrinkage (mash) combinatorial posterior probabilities that a gene-trait association (GTA) is present in class-level cell-type A and absent in cell-type B (). Only gene-trait pairs with \(P({\mathrm{Cell}}_{A}\cap {\mathrm{Cell}}_{B}^{C})\ge 0.6\) in at least one comparison (16 gene-trait combinations) are shown; gray cells indicate combinations without an applicable estimate. Full names of disorders are described in Supplementary Table 12. Full names of all cell-types are defined in Supplementary Table 3. B, Forest plot of BIN1 association effect sizes for Alzheimer’s disease (AD) across cell-types in Million Veteran Program (MVP) individual-level snTWAS (I-snTWAS). The circle represents the effect size, the horizontal bars represent the 95% confidence intervals, and asterisks mark false discovery rate (FDR)-significant associations. Analysis is performed on 37,038 AD cases and 401,544 controls (Supplementary Table 12). C, Top 10 hierarchically pruned pathways enriched for AD snTWAS signals across cell-types, prioritized by Correlation Adjusted MEan RAnk (CAMERA) two-sided p-value; asterisks indicate FDR-adjusted two-sided p-value ≤ 0.05, and at most two significant pathways per cell-type are displayed to increase cell-type-specific representation. D, Analogous pathway enrichment results for schizophrenia (SCZ), showing the top 10 hierarchically pruned pathways across cell-types with the same plotting conventions as in panel C. Statistical test, p-value correction, etc. are performed identically to panel C.

Source data

Extended Data Fig. 4 Relationship of Alzheimer’s disease (AD) summary-level single-nucleus transcriptome-wide association study (S-snTWAS) genes to AD-associated changes in Microglia.

A, Relationship between AD S-snTWAS z-scores and Multiscale Embedded Gene Co-expression Network Analysis (MEGENA) co-expression modules derived from freshly isolated primary human microglia. The top heatmap shows Pearson’s correlations between module enrichment scores for AD S-snTWAS z-scores (x-axis) and module enrichment scores for AD-related neuropathology measures (Clinical Dementia Rating (CDR), Braak stage, and beta-amyloid plaque density; y-axis) across all 306 modules containing at least 50 genes. Asterisks mark false discovery rate (FDR)-adjusted two-sided p-value ≤ 0.05. The bottom heatmap shows, for each cellular population (snBulk and class-level analyses), the number of MEGENA modules significantly enriched (FDR-adjusted two-sided p-value ≤ 0.05) for the corresponding S-snTWAS signature. B, Heatmap for the Correlation Adjusted MEan RAnk (CAMERA) enrichment of neuropathology and S-snTWAS signatures in the 3 FDR significant modules from panel A (M406, M824, M947). Asterisks indicate FDR-adjusted two-sided p-value ≤ 0.05, and white text is used in selected cells to improve legibility. The color corresponds to the -log₁₀(two-sided p-value) of CAMERA analyses multiplied by the sign of the observed enrichment, with negative values reflecting the enrichment of the modules for genes negatively associated with the indicated measures. C, Heatmap of top pathways enriched in the 3 FDR-significant MEGENA modules. Asterisks indicate fisher’s exact test FDR-adjusted two-sided p-value ≤ 0.05. Pathways were renamed for greater clarity of the underlying biology. Upward arrows denote upregulation and downward arrows denote downregulation in the differentially expressed gene signatures derived from the Molecular Signatures Database (MSigDB) C4 Immune gene set collection. Aliases for each pathway are listed in Supplementary Table 24. D, MEGENA co-expression network for module M947 showing gene-gene connectivity; node color encodes AD S-snTWAS significance and node shape denotes MEGENA key driver genes.

Source data

Supplementary information

Supplementary Information (download PDF )

Supplementary Notes 1–7, Supplementary Figs. 1–32, descriptions for Supplementary Tables 1–42 (tables supplied separately), legends for External Data 1–4, Supplementary References and Supplementary Acknowledgements.

Reporting Summary (download PDF )

Supplementary Tables (download ZIP )

Supplementary Data (download ZIP )

Source Data for Supplementary Figs. 1–32.

Peer Review file (download PDF )

Source data

Source Data Figs. 1, 2, 3, 4, 5, 6 and Source Data Extended Data Figs. 3, 4 (download ZIP )

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, 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 changes were made. 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/4.0/.

Reprints and permissions

About this article

Cite this article

Venkatesh, S., Kosoy, R., Wu, Z. et al. Single-nucleus transcriptome-wide association study of human brain disorders. Nature 657, 1016–1026 (2026). https://doi.org/10.1038/s41586-026-10836-6

Download citation