Single-cell atlas of transcriptomic vulnerability across brain disorders

The human brain is a highly complex organ composed of billions of functionally diverse cells. Under pathogenic stress, cellular and molecular responses are often convoluted and contextual, so understanding their dysfunction in disease is challenging. Single-cell approaches that facilitate analysis of the molecular changes that occur within individual cells have been particularly helpful in understanding the interplay between the different cell types found within complex tissues, as well as the roles of those cells in various disease contexts. In AD, a thorough exposition of cellular heterogeneity in the brain has revealed that studying the coordinated interactions between neurons and glia, and the selective depletion of vulnerable inhibitory neuronal subtypes, is critical for understanding AD pathology^1,2.

Building a large-scale disease atlas at the single-cell resolution creates an opportunity to understand molecular responses at the cellular level and estimate population-level variation in the brain transcriptome. A large sample size provides the power needed to establish robust basal-level conditions and to sufficiently capture the full spectrum of disease pathology. By leveraging cross-disease atlases, studies have revealed shared and distinct patterns of gene expression perturbations in major psychiatric diseases, as well as shared genetic factors leading to molecular convergence^3,4. Thus, characterizing shared transcriptomic vulnerabilities and pathophysiology together has considerable implications for early treatment and the development of effective therapeutics.

Here we introduce the PsychAD cohort, which consists of 1,494 unique brain donors, including neurotypical controls and individuals affected by various neurodegenerative and neuropsychiatric diseases (NPDs). Focusing on the dorsolateral prefrontal cortex (DLPFC), we generated a single-nucleus RNA sequencing (snRNA-seq) dataset comprising over 6.3 million nuclei representing 27 distinct subclasses of cells. The dataset is sufficiently powered to identify molecular signatures from multiple traits while accounting for individual variation. To support interpretation and enable future analyses, we evaluated power and effect size variability across cell types (Methods and Supplementary Fig. 5a–g). We identify shared transcriptomic vulnerability, defined as overlapping patterns of gene expression changes across diseases, which point to common pathogenic mechanisms. These shared signatures provide insights into potential molecular mechanisms of early intervention that go beyond traditional disease boundaries. Furthermore, we find that shared transcriptomic vulnerability between disorders is concordant with genetic co-heritability, suggesting that common genetic risk factors partly drive shared gene expression changes. Deep phenotyping of AD trajectories using tau pathology and clinical dementia status suggests a potential link between immune and brain vasculature dysfunctions. In summary, our study provides a rich and comprehensive resource for exploring the cellular and molecular mechanisms of brain function and dysfunction across multiple neurodegenerative diseases and NPDs.

The PsychAD cohort comprises 1,494 unique donors (Fig. 1a, Supplementary Fig. 2 and Supplementary Table 1). Brain tissue specimens were obtained from three sources: 1,042 donors from Mount Sinai NIH Neurobiobank (MSSM), 300 from Human Brain Collection Core (HBCC) and 152 from Rush Alzheimer’s Disease Center (RADC). The cohort covers the whole lifespan of postnatal ages between 0 and 108 years, roughly equal numbers of male (n = 723) and female (n = 771) individuals, and represents a diverse range of disease phenotypes, including AD, DLBD, Vas, tauopathy (Tau), frontotemporal dementia (FTD), PD, schizophrenia (SCZ) and bipolar disorder (BD). The cohort covers a diverse genetic background and over 30% of the donors (n = 509) were of non-European (EUR) ancestry.

We streamlined the unified processing of the data, as well as the harmonization of clinical and technical metadata (Fig. 1b, Supplementary Fig. 1a and Methods). Frozen brain specimens were randomized and processed in batches of six. Equal numbers of nuclei from each sample were pooled together and each pool was processed for two independent snRNA-seq reactions to generate a technical replicate. After quality control (Supplementary Fig. 1c–i), the final datasets consisted of 6,320,459 nuclei. To characterize transcriptomic vulnerability across multiple neuropsychiatric and neurodegenerative diseases (NDDs), we organized the analysis into three tiers.

First, we performed cross-disorder analyses (Fig. 1c). To ensure analytical power and interpretability, we focused on eight well-represented disorders with clear pathological definitions and sufficient donor representation, including balanced numbers of cases and controls. In detail, using a subset of 1,160 donors, at least 17 years of age, with minimal comorbidities, we targeted six NDDs (AD, DLBD, Vas, Tau, FTD and PD) and two NPDs (SCZ and BD). To estimate the sharing of transcriptomic vulnerability, we compared donors affected with NDDs and NPDs against the baseline of 319 neurotypical controls.

The second tier of analysis focused on the stage of AD progression based on pathological and cognitive impairment measures using a subset of 696 individuals (Fig. 1d). To identify cell-type-specific roles in disease onset and disease trajectory, we compared two characteristic neuropathological abnormalities of the AD brain: the accumulation of amyloid-β (Aβ) plaques, measured using the CERAD plaque density score⁵, and tau-based neurofibrillary tangle (NFT) pathology measured using Braak staging^6,7,8, along with cognitive impairment.

In the third tier, we surveyed NPSs within 234 individuals with AD (Fig. 1e and Supplementary Fig. 1b). NPSs are core features of AD and are common in patients with dementia⁹. We broadly categorized NPSs into three groups—(1) depression or mood-related, (2) weight loss or psychomotor agitation (PMA), and (3) weight gain, insomnia or suicidal ideation—on the basis of co-occurrence estimates, sharing molecular mechanisms that lead to increased prevalence with disease severity.

To understand heterogeneous human cortical tissues in disease contexts, we require a cell type taxonomy that is robust to ageing, disease phenotypes and various sampling and technical biases. After unified computational processing, quality control and batch normalization of snRNA-seq libraries representing 1,494 dissections, processed in duplicate (Methods), we annotated cell types of the human DLPFC using the cell taxonomy of the primate DLPFC¹⁰ and human primary motor cortex¹¹ as baseline references. The resulting taxonomy was organized into three hierarchical levels, identifying eight broad cell classes, 27 subclasses and 65 functionally distinct subtypes (Fig. 2a and Supplementary Table 2). Each level of the annotation hierarchy represents a slice in the clustering dendrogram. At the top, the class level defines eight major cell types, including two broad neuronal cell types: glutamatergic excitatory (EN) and GABAergic inhibitory neurons (IN); three glial: astrocytes (Astro), oligodendrocytes (Oligo) and oligodendrocyte progenitor cells (OPCs); and three non-neuronal cell types: immune cells (Immune), mural and vascular cells (Mural) and endothelial cells (Endo). Subsequent levels, namely subclass and subtype, were derived by iteratively reclustering the subset of cells by gene matrix using a new set of variable genes relevant to the particular cell type (see the ‘Defining cellular taxonomy using iterative clustering’ section of the Methods). The cellular taxonomy was consistent and well represented across all three brain sources (Supplementary Fig. 3a,b). Our most granular-level subtype annotation was robust and invariant to donor and technical variables (Supplementary Fig. 3c). Neuronal cells made up 38.4% (EN 23.0% and IN 15.4%) of all cells, with oligodendrocytes being the next most abundant, at 36.1%. Major cell types and subtypes matched well when compared to previous DLPFC cellular taxonomies^1,10 (Supplementary Fig. 3d–f and Supplementary Table 11).

The subclass level distinguished the EN and IN classes into 10 and 7 subclasses, respectively (Supplementary Fig. 3g,h). Different types of neurons, especially ENs, are organized into six horizontal layers (L1–L6) that are distinct in both cytoarchitecture and function¹². In situ spatial transcriptomics data were used to confirm that the EN subclasses were spatially distinct and found in their respective neocortical layers (Fig. 2b,c and Supplementary Fig. 4j). The EN subclasses were denoted by their laminar organization (L2–6) and axon projection characteristics (intratelencephalic (IT), extratelencephalic (ET), near projecting (NP), corticothalamic (CT) and L6B). IN subclasses were determined using their characteristic marker genes (Fig. 2d and Supplementary Fig. 3i); LAMP5⁺LHX6⁺ ivy cells (IN_LAMP5_LHX6), LAMP5⁺RELN⁺ neurogliaform cells (IN_LAMP5_RELN), parvalbumin⁺ basket cells (IN_PVALB), parvalbumin⁺ chandelier cells (IN_PVALB_CHC), somatostatin⁺ Martinotti and non-Martinotti cells (IN_SST), VIP⁺ inhibitory neurons (IN_VIP) and ADARB2⁺ homologues of mouse Sncg⁺ inhibitory neurons¹⁰ (IN_ADARB2). Consistent with previous reports^10,13,14, IN subclasses did not display the same degree of laminar organization as EN subclasses (Supplementary Fig. 4i). A more careful examination of the spatial density of nuclei from each subclass within coarsely defined laminar regions (Methods) revealed that IN distribution is not uniform across cortical layers. For example, we observed a preference for the ADARB2⁺ and VIP⁺ interneurons for outer cortical layers, while both PVALB⁺ interneuron subclasses were distributed across layers L2–6 with less prominent bias (Supplementary Fig. 4e,f). Notably, we observed SST⁺ inhibitory neuron subtypes had distinct layer-specific enrichment patterns (Supplementary Fig. 4f,k). Some previously annotated rare inhibitory neuron types, such as SST⁺NPY⁺ or SST⁺HGF⁺, were not distinguished at the subclass level, but were identified at the subtype level.

As expected, neuronal subclasses exhibit distinct functional characteristics from those of non-neurons (Fig. 2e). ENs are enriched in pathways associated with synaptic vesicle priming and neurotransmitter secretion, whereas INs are enriched with functions related to receptor signalling and ion transport. Immune cell types show unique enrichments in functions related to cytotoxic immune responses, while mural cells are involved in vessel morphogenesis. It has been shown that cell types are differentially implicated in mediating disease risk^15,16. We therefore further annotated cell types with disease genome-wide association studies (GWAS) using single-cell disease-relevance score (scDRS) and found that neurological diseases, including NDDs such as AD and neuroimmune diseases such as MS, largely involve immune and glial cell types, whereas psychiatric diseases are predominantly associated with neurons (Fig. 2f and Supplementary Fig. 3j).

Transcriptomic variation is influenced by genetic differences among individuals, phenotypes and technical factors. To understand the determinants of the overall transcriptome variation in our cohort, we used population-scale snRNA-seq data to partition gene expression variance by cellular variables (including cell type and fraction of mitochondrial and ribosomal genes), donor-level variables (including individual ID, age, sex, genetic ancestry and diagnosis) and technical variables (including source of the tissue specimen, post-mortem interval, technical replicates and sequencing depth) (Fig. 3a and Methods). Across all genes, a mean of 50.5% of the total expression variance can be attributed to cell type variation, while interindividual variation explains 7.5% and unexplained residual variation accounts for 41.1%. The remaining variables explained, on average, less than 1% of the total variance. We prioritized genes with the highest expression variation across variables (Fig. 3b and Supplementary Fig. 6b–f). Among the genes with high cell type variation (including ATP8A2, DOCK3 and GABRB3), most are differentially expressed between neuronal and non-neuronal cell types (Fig. 3c). As expected, genes varying across sexes were located in sex chromosomes (Supplementary Fig. 6d), while those with high variation across tissue sources were mostly mitochondrial genes (Supplementary Fig. 6e), possibly due to technical differences in physiological and environmental factors from dissection and handling at the respective tissue sources^17,18. Genes with high variation across diagnosis were implicated in stress-response and housekeeping functions (Supplementary Fig. 6a,f). The top gene, cold-inducible RNA-binding protein (CIRBP), responds to hypothermic stress and protects cells from hypoxia-induced neurotoxicity and its expression has been shown to be inversely associated with patient survival in cancer^19,20,21. Note that, in addition to its variation across diagnoses, CIRBP also had the highest variation with PMI. As this approach examines expression variation shared across cell types, the low fraction of variance explained by diagnosis indicates the need for cell-type-specific analysis.

Notably, we observed that interindividual differences explained 82.2% of the variation in expression of the ADP ribosylation factor-like GTPase 17B (ARL17B) gene (Fig. 3d). Several adjacent genes, such as KAT8 regulatory NSL complex subunit 1 (KANSL1) and ARL17A (a paralogue of ARL17B), also had high interindividual variation (Supplementary Fig. 6c). ARL17B and KANSL1 are localized in the disease-associated MAPT locus (17q21.31), are often found as a fusion transcript (KANSL1::ARL17B), frequently undergo polymorphic translocation^22,23 and have been implicated in neurological disorders, including ALS, PD and MS^24,25,26,27. To examine possible genetic causes for variation in expression, we examined normalized gene expression at the donor level, stratified by MAPT haplotypes (Fig. 3e and Methods). We observed two distinct patterns of ARL17B expression that could be potentially linked to H1 and H2 MAPT haplotypes, with lower expression associated with the H1H1 genotype. We observed stratification of haplotypes by genetic ancestry, where the H2H2 genotype was almost exclusively found within EUR ancestry, consistent with previous reports^28,29 (Supplementary Fig. 6g). We also replicated a previous finding indicating that the H1 haplotype is associated with PD susceptibility, and observed that the H1H1 genotype increases PD risk with an odds ratio (OR) of 4.125 (Supplementary Fig. 6h; P ≤ 0.0273), albeit much higher than previously reported (1.42 (ref. ³⁰) or 1.46 (ref. ³¹). Moreover, we tested the contribution of the H1 haplotype to AD among non-APOE4 carriers³² but did not find a significant association (P ≤ 0.302). The H1 haplotype exhibited cis-regulatory effects on several genes near the MAPT locus that are relevant to neurodegeneration³³ (Supplementary Fig. 6i). Specifically, H1 haplotypes in neurons and astrocytes were associated with increased expression of leucine rich repeat containing 37A genes (LRRC37A2/A3). LRRC37A proteins have been shown to interact with α-synuclein and co-localize with Lewy bodies in PD brains³⁴. Another gene upregulated by the H1 haplotype was pleckstrin homology domain-containing family M member 1 (PLEKHM1), which has been linked to PD susceptibility by GWAS^35,36. The H1 haplotype also demonstrated interchromosomal effects, notably on phosphodiesterase 8B (PDE8B), a gene implicated in movement disorders^37,38,39,40. Furthermore, heterogeneous nuclear ribonucleoprotein K (HNRNPK), a gene involved in synaptic plasticity regulation⁴¹, was downregulated by the H1 haplotype.

Genes with lower interindividual variation had higher genetic constraints, as measured by the gnomAD⁴² loss-of-function observed/expected upper-bound fraction (LOEUF) score or GeneBayes⁴³ constraint metric (S_het) (Fig. 3f and Supplementary Fig. 6l). These conserved genes were functionally implicated in nervous system development and synaptic functions, indicating that the gene program in neuronal development is conserved across individuals (Supplementary Fig. 6m). On the other hand, genes with high interindividual variability were enriched with basal cellular functions (that is, housekeeping genes, RNA metabolism and, in particular, cellular response to heat stress), corroborating previous findings^44,45 (Fig. 3g and Supplementary Fig. 6a,m).

Using the 318 neurotypical donors as a baseline, we systematically evaluated variation in the cellular composition of the DLPFC across eight different NDDs (AD, DLBD, Vas, Tau, PD and FTD) and NPDs (SCZ and BD). Using all subclass-level cell types, we found that the overall cell type composition changes were broadly stratified by NDDs and NPDs, and that they formed distinct clusters (Fig. 4a). On average, correlations were greater when focusing solely on neurons, underscoring their critical role in the aetiology of neurological diseases. Notably, stronger similarities were observed between DLBD and both AD and Vas. Moreover, similarities between FTD–AD (all cells) and FTD–Tau (neurons) were observed. Exploring each subclass further, we identified a notable overlap in the prevalence of neuronal and glial cell types within the same class of diseases (Fig. 4b and Supplementary Table 3). Specifically, NDDs were characterized by a higher abundance of non-neuronal cells, particularly VLMCs and other vascular cell types (Supplementary Fig. 7c,e,f), as well as elevation of specific inhibitory neuron subclasses, namely IN_LAMP5_RELN, IN_LAMP5_LHX6 and IN_ADARB2. By contrast, NPDs were predominantly associated with an increase in neuronal cells, particularly deep-layer ENs in L5–6. To further identify specific subtypes responsible for driving the compositional changes in vascular cell types, we used subtype-level annotation to analyse compositional variation in 8 NDDs and NPDs (Fig. 4c and Supplementary Fig. 7d). From this, we identified dominant subtypes of each subclass that further differentiated NDDs and NPDs. For example, vascular leptomeningeal cells (VLMCs) are barrier-forming fibroblasts of the brain⁴⁶, and they are transcriptionally segregated into three subtypes; two meningeal VLMCs (VLMC_DCDC2 and VLMC_SLC4A4) and one perivascular VLMC (VLMC_ABCA6). Our subtype-level analysis indicates a polarized response of meningeal VLMC_DCDC2, where their increased proportions are specifically associated with most NDDs.

Analysing gene expression signatures across multiple disorders can reveal shared molecular pathways, facilitating improved diagnostics and targeted therapeutics. Thus, we sought to assess the extent of sharing across NDDs and NPDs. Using a pseudobulk approach and precision-weighted regression modelling for large-scale single cell datasets, implemented in the dreamlet software⁴⁷ (Methods), we performed a comprehensive analysis of differentially expressed genes (DEGs) across eight different NDDs and NPDs. After identifying DEGs for each disease contrast and tissue source, we performed a standard frequentist fixed-effect meta-analysis across the three tissue sources to define consensus disease trait signatures. We then applied a multivariate Bayesian meta-analysis^47,48 across cell types and disease traits to produce posterior estimates of the effect size and the posterior probability that each effect is non-zero. Using a composite testing approach, we used these results to estimate the posterior probabilities that each effect is specific to a trait or shared across traits. We then evaluated, in a cell-type-specific manner, cross-disease sharing and decomposed the total disease effects into shared and disease-specific components (Supplementary Table 4). We demonstrate that cross-disorder signatures have a higher degree of disease-effect sharing (Supplementary Fig. 8a) and encompass crucial transcriptional processes, such as mRNA splicing and processing, as well as protein localization to mitochondria (Fig. 5a and Supplementary Fig. 8b). The observed shared signatures affect genes that are critical for proper cell function align with the omnigenic model of disease risk⁴⁹, and further support the pleiotropy of these genes, which influence multiple disorders both genetically and transcriptionally. Genes with disease-specific effects exhibited a higher burden of rare variants compared with those with shared effects (Supplementary Fig. 8e), with statistically significant enrichments observed in VIP⁺ and PVALB⁺ interneurons. While this supports the relevance of disease-specific components to core disease biology⁵⁰, the signal was limited to a subset of cell types. To facilitate the identification of core disease-relevant functional genes, we discounted genes with shared effects from the overall DEG expression profiles. Using the residual disease-specific effects, we quantified the pairwise transcriptomic similarity between traits (Fig. 5b). Similarities between pairs of NDDs or NPDs were greater compared with the NDD–NPD contrast, with AD, DLBD, Vas and PD being the most similar. Subsequent meta-analysis using disease-specific effects of these four transcriptomically similar traits implicated neuronal development and synaptic signalling pathways in interneurons (IN_LAMP5_RELN, IN_ADARB2 and IN_PVALB), as well as vasculature development in the VLMC subclass (Fig. 5c and Supplementary Fig. 8c).

Comparison of pairwise transcriptome similarity against pairwise trait co-heritability (common genetic variation) can estimate the contribution of shared genetic influences on molecular phenotypes. We calculated the shared heritability (r_g) between disease pairs using GWAS and found that pairs within NDDs or NPDs exhibited a higher degree of shared heritability (that is, DLBD–PD, SCZ–BD) compared with pairs spanning NDDs and NPDs (Fig. 5d). Quantifying the pairwise transcriptome similarity (ρ_s) across NDDs and NPDs revealed varying degrees of transcriptomic overlap, with some exhibiting high similarity, particularly in disease pairs such as AD and DLBD (Fig. 5e). By correlating the shared heritability against the average transcriptome similarity across cell types, we obtained a heritability–transcriptome concordance (r_g − ρ_s), which indicated a strong positive correlation (Spearman’s ρ = 0.604) between genetic and transcriptomic similarities (Fig. 5f and Supplementary Fig. 8f). The results show that diseases with higher genetic overlap also tend to have more similar transcriptomic dysregulation, suggesting that shared transcriptomic vulnerabilities are coupled with underlying genetic architecture, an unsurprising result that corroborates previous observations^3,51,52. We next extended the r_g − ρ_s comparison by considering the transcriptional concordance for each cell type (Fig. 5g and Supplementary Fig. 8d). Among all cell types, Chandelier cells (IN_PVALB_CHC), a subset of GABAergic interneurons, had the highest concordance with the pairwise trait heritability. In summary, our approach enabled us to dissect disease signatures into shared and distinct components, revealing that core disease signatures have significant overlap in gene expression and genetic risk across NDDs and NPDs. These findings further our understanding of disease mechanisms in the brain, paving the way for novel therapeutic strategies targeting shared or distinct pathways.

Next, we focused on characterizing the transcriptomic variation in AD, including analysis of case–control comparisons and different phenotypes that capture disease severity, such as plaque density (that is, CERAD scores), NFT progression (that is, Braak stage) and level of cognitive impairment (such as CDR). The variation in cell type composition reveals distinct patterns associated with AD pathologies compared with normal ageing (Fig. 6a, Supplementary Fig. 9a and Supplementary Table 5) and was consistent with previous studies^1,2 (Supplementary Fig. 10a,b). Furthermore, we replicated our findings using a published MERFISH dataset², showing that the observed change is robust to the technology used (Supplementary Fig. 10c). Comparing compositional shifts during normal ageing and AD (Fig. 6b), we identified an AD-specific increase in smooth muscle cells (SMCs), a vascular cell type that normally declines with age, suggesting an AD-specific vulnerability in these cells. Likewise, changes in most excitatory neurons were discordant, suggesting that the mechanisms that led to neuronal loss were AD specific. Notably, changes more pronounced in AD, such as the increase in LAMP5⁺LHX6⁺ inhibitory neurons (IN_LAMP5_LHX6) and the loss of SST⁺ inhibitory neurons (IN_SST) with higher CERAD scores indicates a potential association between this neuronal subtype and Aβ pathology. Moreover, the loss of L2–3 IT neurons in AD (based on case–control and CERAD comparisons), which is not evident in normal ageing, suggests a specific vulnerability for this neuronal subtype to AD. We next experimentally validated the compositional shifts using RNAscope (Methods), confirming the decline in L2–3 IT excitatory and SST⁺ inhibitory neurons and the prevalence of VLMCs in AD (Supplementary Fig. 9e,f). Taken together, our findings highlight the cell type specificity of AD vulnerability.

To better understand transcriptomic variation in AD, and to identify the cell subclasses most affected, we analysed the changes in cases of AD due to co-occurring conditions. It is estimated that more than 80% of patients with AD will exhibit at least one NPS over the course of their illness that substantially impacts their clinical outcome^9,53, suggesting that there are at least some shared molecular mechanisms between serious mental illness and AD. To evaluate the cell types associated with NPS prevalence, we applied compositional variation analysis to three categories of NPS based on co-occurrence estimates. Notably, using age-matched groups of patients with AD with or without NPSs (Supplementary Fig. 9b), we found that patients with AD experiencing weight loss and PMA have an increased ratio of excitatory neurons, especially deep-layer neurons in L5–6 (Fig. 6a and Supplementary Fig. 9a,c,d). Our results corroborate previous findings⁵⁴ that changes in deeper PFC layers (L5–6) are associated with certain types of NPS.

After identification of vulnerable cell subclasses in AD, we questioned whether their roles in AD pathology are damaging, protective, causal or derived. As demonstrated in previous studies⁵⁵, we used mediation analysis to decipher the causal relationships between various cascades of events leading to AD onset and progression. We performed the analysis on the base hypothesis that polygenic risk (AD PRS; Methods) contributes to plaque accumulation (mean density of neuritic plaques), which in turn affects tau progression (Braak), leading to dementia. We tested 14 subclasses of which the cell type composition significantly changed in any of the AD contrasts (false-discovery rate (FDR) < 0.05) and identified that microglia, VLMCs, SST⁺ and LAMP5⁺LHX6⁺ inhibitory neurons exhibit significant average causal mediation effects (ACME, P < 0.05) (Fig. 6c and Supplementary Table 13). Beyond the direct effects of polygenic risk on tau pathology, we observed significant indirect effects mediated by plaque accumulation (ACME = 0.0903, P < 2 × 10⁻¹⁶). Tau progression results in more VLMCs mediated by an increase in microglia (ACME = 0.00477, P = 0.005). Conversely, an increase in microglia leads to a decrease in LAMP5⁺LHX6⁺ inhibitory neurons. Such changes in the levels of both VLMCs (ACME = 0.00515, P = 0.0002) and LAMP5⁺LHX6⁺ inhibitory neurons (ACME = −0.00269, P = 0.031) appear to exacerbate dementia. Furthermore, plaque accumulation is mitigated by SST⁺ inhibitory neurons (ACME = 0.1475, P < 2 × 10⁻¹⁶). Increased microglia lower SST⁺ inhibitory neurons (coef = −0.18, P = 5.04 × 10⁻⁷), as does increased tau (coef = −0.06, P = 1.82 × 10⁻⁴), contributing to more plaque accumulation (coef = −2.01, P = 7.47 × 10⁻⁶).

To determine how gene programs change in response to increasing severity of AD pathology, we characterized the DEGs in 27 subclasses (Fig. 6d, Supplementary Fig. 11b and Supplementary Table 6). In general, DEGs had high concordance across different AD pathology variables, consistent with previous reports¹, and the disease signatures were highly concordant with previous studies^1,2 (Supplementary Fig. 12a–e). Overall, DEGs in AD (FDR < 0.05) can be summarized as upregulation of genes in vascular cell classes (Mural and Endo), while genes in neurons tend to be downregulated (Supplementary Fig. 11a). We observed that microglia-specific DEGs, including DPYD, IL15 and PTPRG^47,56,57, generally exhibited higher effect sizes compared with those of other subclasses. Further characterizing microglia-associated gene signatures, we found that they are specifically enriched for pathways involved in negative regulation of cell motility and migration, as well as response to lipoproteins (Fig. 6e and Supplementary Fig. 11c). Gene expression changes in neurons were largely related to synaptic functions, including synapse assembly, development, signalling and membrane transport. Lastly, genes altered in VLMCs were implicated in muscle tissue development.

Recent studies highlight the complex, time-dependent changes in gene expression and cellular function during AD progression. Early stages are characterized by increased inflammation, astrocyte reactivity, blood–brain barrier dysfunction and heightened vulnerability in select neuronal populations^{1,2,46,55,58,59}. These are followed by a more dysregulated immune response and broader neuronal loss in the late stages^{60,61,62,63,64,65}. To further investigate the gene expression dynamics underlying AD pathogenesis, we used neural network models to generate two independent cell-type-specific disease trajectories from semi-quantitative measures of AD progression (that is, tau proteinopathy and cognitive decline) (Fig. 7a). As tau proteinopathy and dementia are correlated, models that do not disentangle their effects may miss genes with opposing expression patterns. For example, if a gene’s expression rises with tau but falls with dementia, these effects could cancel one another out, obscuring significant changes. Thus, to decorrelate Braak and dementia model predictions, we equally sampled all combinations of Braak stage and dementia status during model training (Methods and Supplementary Fig. 14; potential consequences of not decorrelating are shown in Supplementary Fig. 15). The accuracies of the Braak and dementia model predictions were significantly above chance for all eight cell classes (P < 1 × 10⁻⁴, bootstrap) (Fig. 7b), with predictions generated on cross-validated data, using non-overlapping donors for training and testing. Visualizing the model predictions of Braak (hereafter, disease pseudotime) in a UMAP space for microglia shows that actual Braak values closely align with the model predictions (Fig. 7c and Methods).

Many past studies have suggested that gene expression can evolve nonlinearly with disease progression⁶⁶. Thus, we performed principal component analysis (PCA) analysis of the concatenation of all disease pseudotime trajectories (Fig. 7a) across all genes and cell types, revealing the presence of a nonlinear transition between early and late disease stages (Supplementary Fig. 16), and showing that glial and vascular trajectories were more nonlinear than neuronal trajectories (Fig. 7d and Supplementary Fig. 17; P < 0.05, Wilcoxon rank-sum test). In addition to the disease pseudotime axis, we also calculated a dementia resilience score that measured how gene expression is correlated with predicted dementia, conditioned on the predicted disease pseudotime (Methods). A gene is considered protective if increased expression is associated with a decrease in predicted dementia, conditioned on the disease pseudotime, and damaging if increased expression is associated with an increase in predicted dementia. We then determined the top 32 Gene Ontology (GO) Biological Process (BP) pathways that best summarize the biological changes associated with early and late changes in disease pseudotime and resilience (Fig. 7e, Supplementary Figs. 18–20 and Supplementary Table 7).

A semantic clustering of the pathways identified five functional clusters. The first cluster was primarily related to synaptic function, characterized by downregulation in OPCs in the early stages of AD, and then upregulation in neurons and downregulation in immune cells in later stages. For these cell types, increased expression was associated with diminished resilience (that is, damaging). A visualization of these damaging synaptic pathways is shown in Supplementary Fig. 21a,b. The late increase in genes implicated in these pathways could be linked to compensatory mechanisms after synaptic loss⁶⁷. The second cluster was related to cell metabolism pathways implicated in protein translation^2,68, mitochondrial function^2,69 and acidification⁷⁰, which showed strong downregulation in neurons with increasing disease pseudotime, consistent with past studies⁷¹. Decreased expression in these metabolic pathways was associated with decreased resilience (increase in metabolism was protective). The third cluster was related to cell stress, including pathways related to chaperone-mediated protein assembly, autophagy and response to heat. In non-neuronal cells, both early and late increases in cell stress were strongly associated with cognitive decline. The fourth cluster was related to immune response and inflammation. Immune cells were implicated immune responses such as antigen presentation and helper T cell response, as well as monocyte differentiation, all of which were associated with dementia resilience. A visualization of these protective immune pathways is shown in Supplementary Fig. 21c,d. Lastly, vascular cells were implicated in damaging changes associated with adherens junction assembly, response to interleukin-17 (IL-17), glycolysis, endothelial cell differentiation and nitric oxide synthesis.

While this analysis included all pathways, we also aimed to identify potentially causal pathways driving disease pathogenesis. To do so, we repeated our analysis, focusing only on GO BP pathways enriched in risk genes (Methods). One particularly notable example is the negative regulation of peptidyl-threonine phosphorylation (Supplementary Figs. 22 and 23), which has been implicated in tau phosphorylation⁷². The mean expression of this pathway selectively decreases early in both neural cell types and includes the risk gene SPRED2 (ref. ⁷³).

As previous studies have highlighted the importance of the immune response in AD^{57,61,62,63,64,74,75,76}, we examined this cell class in more detail (Fig. 7f–i). To help visualize the nonlinear dynamics of the immune response in AD, we categorized gene expression trajectories into four modules on the basis of their response to increasing tau proteinopathy: early increasing, early decreasing, late increasing and late decreasing (Fig. 7f). Gene enrichment of these four trajectory modules showed that pathways involved in macrophage colony stimulation and pyrimidine metabolism were all upregulated in the early stage (Fig. 7e), whereas migration, purinergic signalling and integrin activation were downregulated. In the later stages of disease, lipid-related pathways such as adipose tissue development and triglyceride metabolism, in addition to cytokine production, increased, while synapse-related pathways were downregulated. Gene-set enrichment analysis using GWAS summary statistics for the top 250 genes in each trajectory module revealed that late-increasing genes were the most strongly associated with AD (FDR ≤ 1.7 × 10⁻⁷), although early-increasing (FDR ≤ 1.7 × 10⁻³) and early-decreasing (FDR ≤ 7.9 × 10⁻⁴) genes were also significantly associated (Fig. 7h and Supplementary Table 8). Notably, the late-increasing trajectory module was also significantly correlated with other NDDs, including MS and PD. We also find that changes across many cell classes significantly overlap with NPD risk loci (Supplementary Fig. 24 and Supplementary Table 8).

We further examined the genes associated with trajectory modules (Fig. 7i; the top genes for all cell classes are listed in Supplementary Table 9). Several markers for homeostatic microglia⁵⁷, including CX3CR1, NAV2 and P2RY12, were among the top 5 early decreasing genes with increasing disease pseudotime (FRMD4A and CERC2 were ranked in the top 20 out of 17,265 coding genes). By contrast, several early-increasing genes, such as ACSL1, DPYD and CD163 were recently implicated in a pathogenic lipid-droplet accumulation phenotype in patients with APOE4/4 AD⁶¹ (another implicated gene, NAMPT, ranked 134th; Supplementary Fig. 26a and Supplementary Table 9).

To further resolve the immune response at a finer time resolution, we performed pathway enrichment using a sliding window across disease pseudotime (Supplementary Fig. 25). Our results suggest two distinct immune phases⁶⁰: an initial innate response, followed by a more adaptive and monocyte-driven response, along with increased cytokine production. Lipid homeostasis and efflux pathways, many involving known AD risk genes, increased throughout disease progression, particularly in later stages (Supplementary Figs. 25 and 26b), supporting the hypothesis that microglia develop a lipid-droplet accumulation phenotype that may exacerbate disease progression^61,62,63.

Here we report a comprehensive single-nucleus transcriptomic atlas of the human DLPFC, generated using 1,494 donors, including neurotypical controls and individuals affected by a range of complex neurodegenerative and/or neuropsychiatric conditions. Our analyses provide insights into the cellular heterogeneity of the DLPFC, including the observation that about 10% of total transcriptomic variation can be attributed to interindividual differences. Notably, we find that genes with lower interindividual variability are under higher genetic constraints. These gene programs are enriched with neuronal development pathways, indicative of their critical roles in development and fitness.

We used our disease atlas to characterize cross-disease and disease-specific responses to pathologic conditions, revealing shared and distinct cellular composition profiles among NDDs and NPDs. Shared signatures include genes critical for proper cellular function (mRNA processing and protein localization), and align with the omnigenic model⁴⁹, in which a portion of heritability in complex traits can be explained by effects on peripheral, often regulatory, genes that can have indirect, subtle and cumulative roles in disease. By focusing on core disease-specific functional genes, we found that disease pairs with higher genetic risk overlap tend to have greater cell-type-specific transcriptomic concordance, corroborating previous observations from bulk assays³ that genetic factors contributing to disease susceptibility can influence transcriptomic alterations in similar ways. The concordance between shared differential expression patterns and genetic co-heritability suggests that overlapping transcriptomic signatures reflect core biological processes influenced by pleiotropic risk variants. This concordance provides insight into molecular pathways that may be causally linked to disease aetiology across disorders. Cell types with higher concordance, including deep-layer excitatory neurons and PVALB⁺ Chandelier cells, probably represent critical sites of genetic risk convergence and may warrant prioritization in mechanistic studies. While disease-specific differential expression signatures remain important for understanding disorder-specific pathophysiology, shared signatures offer a window into fundamental processes contributing to comorbidity and pleiotropy.

Among our findings is the observation that the brain vascular system is intricately linked to immune dysfunction in AD. In general, we saw a relative increase in vascular cell types in most NDDs (Supplementary Fig. 6e,f) and demonstrate that levels of VLMCs are particularly elevated in AD among individuals with more severe cognitive impairment (Fig. 6a,c and Supplementary Fig. 9e,f). It has previously been shown that the meningeal lymphatic system has multiple roles in the brain, including waste removal^77,78,79 and the adaptive immune response^46,80. Given the role of VLMCs in these processes^81,82,83,84, their contribution to disease warrants further investigation.

Our analysis of the pathological trajectory of AD aligns with various proposed hypotheses^{71,85,86,87,88,89}, particularly for pathways affected in the earliest stages of the disease. Neuronal, immune and vascular cells exhibit distinct vulnerabilities, with shared alterations across correlated pathways. In neurons, we observed a decrease in metabolic functions (cyclic nucleotide catabolic process, cytoplasmic translation and ATP synthesis) that have been closely linked to synaptic dysfunction and cognitive decline⁷¹. On the other hand, the immune response appears to be generally protective, whereby an early innate immune activation⁸⁵ followed by an adaptive immune response^64,75,76 was associated with resilience to dementia. In particular, monocyte differentiation appeared to have a protective role during both early and late disease stages, consistent with the idea that monocyte-derived macrophages serve as reinforcements when the immune system is under sustained stress⁹⁰. The one exception was the response of immune, mural and endothelial cell classes to IL-17, which was damaging in the early stages of AD. IL-17 has been associated with cognitive decline^91,92 and disruption of the blood–brain barrier⁹³, and blocking IL-17 function has been shown to mitigate cognitive defects in murine models of AD^91,92. Moreover, chaperones have been closely associated with a number of disease related processes, including tau misfolding and aggregation, and neurotoxicity^{86,87,88,89,94,95}; however, their precise roles in AD pathogenesis remains unclear. Our results suggest that chaperones can have opposing roles depending on the cellular context—damaging for glia, immune and endothelial cells but possibly protective for neurons. Lastly, upregulation of endothelial cell stress, differentiation and adhesion/junction pathways were associated with cognitive decline, consistent with the idea that blood–brain barrier dysfunction⁵⁸—a known hallmark of dementia—may trigger compensatory increases in these pathways.

Taken together, these insights deepen our understanding of AD pathogenesis and implicate cellular responses that warrant further investigation. While the human brain comprises numerous regions with distinct functions, in this study, we limited our analysis to the DLPFC. Future studies that include additional brain regions will provide insights into a broader range of phenotypes and offer a more holistic view of the molecular foundations of brain function in both health and disease. Moreover, more detailed spatial analyses, including whole-cell transcriptomics, could further elucidate the physiological contexts that contribute to disease. Analysing larger numbers of cells may also uncover rarer disease-associated cell types or subtypes, while integrating additional omics modalities (for example, proteomics or epigenomics) will provide deeper molecular insights. Overall, the PsychAD single-cell disease atlas serves as a unique and foundational resource to further our understanding of population-level disease-associated transcriptomic variation in the human brain.

Collection and harmonization of clinical, pathological and demographic metadata

Brain tissue specimens were sourced from two brain banks: the Mount Sinai NIH Neurobiobank (MSSM) (1,042 samples) and the NIMH-IRP Human Brain Collection Core (HBCC) (300 samples). Furthermore, samples were obtained from five prospective cohort studies conducted at the Rush Alzheimer’s Disease Center (RADC) (152 samples)^124,125. Thus, the availability of clinical metadata varied as a function of source. We used the following scheme to harmonize available clinical, pathological and demographic metadata: the CERAD scoring scheme for neuritic plaque density⁵ was harmonized for consistency across multiple brain banks, for which the scores range from 1 to 4, with increasing CERAD number corresponding to an increase in AD burden; 1, no neuritic plaque (normal brain); 2, sparse (possible AD); 3, moderate (probable AD); 4, frequent (definite AD). Samples from RADC used consensus summary diagnosis of no cognitive impairment (NCI), mild cognitive impairment (MCI), and dementia and its principal cause, Alzheimer’s dementia^126,127,128. MSSM samples used clinical dementia rating (CDR), which was based on a scale of 0–5; 0, no dementia; 0.5, questionable dementia (very mild); 1, mild dementia; 2, moderate dementia; 3, severe dementia; 4, profound dementia; 5, terminal dementia. After consulting with clinicians, we created a harmonized ordinal variable in which dementia is categorized into three levels of cognitive decline, independently of AD diagnosis: 0, no cognitive impairment; 0.5, MCI; and 1–5, dementia. In addition to AD phenotype, we collected comprehensive demographic (age, sex and genetic ancestry) and technical variables (tissue source, technician, sample batch, postmortem interval (PMI; measured in min), APOE genotype) to describe each cohort (Supplementary Table 1). We described the process for assigning genetic ancestry in a previous data descriptor paper¹²⁹. In brief, we leveraged quadratic discriminant analysis (QDA) to infer genetic ancestry by training our model using data from the 1000 Genomes Project. We used tenfold stratified cross validation to optimize the regularization parameter within QDA¹³⁰, as well as forward selection to identify the optimal number of PCs for genetic ancestry assignments. For samples without genotypic data, we used race/ethnicity as a proxy for inferring genetic ancestry. We emphasize that, while genetic ancestry is a distinct concept from the social constructs of race and ethnicity¹³¹, we leveraged the correlated race/ethnicity variables as proxies to retain those samples in the analyses. Values for superpopulations included EAS, SAS, AFR, AMR, EUR and EAS_SAS, where the category EAS_SAS was assigned to samples with unavailable genotypes with an ‘Asian’ value for race/ethnicity, which can potentially correspond to both EAS and SAS.

Clinical diagnosis of AD

For analysis comparing donors with AD cases and neurotypical controls, a binary clinical diagnosis variable for AD, dx_AD, was defined, as follows: individuals with CERAD 2, 3 or 4, Braak ≥ 3 and CDR ≥ 1 for MSSM or Alzheimer’s dementia for RADC were classified as AD cases. Controls were defined as individuals in the controls_neuropathological_clinical category, where CERAD = {1}, Braak = {0,1,2,3} and secondary diagnosis (including dementia) is not allowed except for MCI.

Measuring AD neuropathology

For analysis comparing donors with pathologic AD, the following variables were used to measure the severity of AD neuropathology: CERAD score⁵. A quantitative measure of Aβ plaque density in which 1 is normal, 2 is possible AD, 3 is probable AD and 4 is definite AD¹²⁶. Braak AD-staging score^6,7,8 measuring progression of NFT neuropathology (Braak and Braak score, or BBScore). A quantitative measure of the regional patterns of NFT density across the brain, in which 0 is normal and asymptomatic, 1–2 indicates initial stages where NFTs begin to appear in the locus coeruleus and the transentorhinal region, 3–4 indicates progression to limbic regions, such as the hippocampus and amygdala, and 5–6 indicates that NFTs are widespread, affecting multiple cortical regions.

Measuring cognitive impairment

For analysis comparing donors with AD-related dementia, the following variable was used to measure the severity of cognitive impairment: clinical assessment of dementia. A harmonized variable of cognitive status based on the CDR scale for MSSM or NCI, MCI and Alzheimer’s dementia for RADC. We used the three-level ordinal categories of clinical dementia to measure the severity of dementia, in which 0 indicates no dementia, 0.5 indicates minor cognitive impairment, and 1.0 indicates definite clinical dementia.

Definition of cross-disorder contrasts

For cross-disorder contrasts, we limited the analysis to any individual with age ≥ 17. Neurotypical controls are defined as any individual CERAD = {1} and Braak = {0,1,2,3} and secondary diagnosis is not allowed. AD is any individual with CERAD = {2,3,4}, Braak = {3,4,5,6} and clinically diagnosed as dementia, and secondary diagnosis not allowed. SCZ is any individual with SCZ diagnosis (SCZ|Schizoaffective_bipolar|Schizoaffective_depressive) and secondary diagnosis not allowed, except for metabolic and eating disorders. DLBD is any individual with DLBD diagnosis (DLBD), and secondary diagnosis can be only AD. Vascular is any individual with vascular diagnosis (vascular) and secondary diagnosis can be only AD. BD is any individual with BD diagnosis (BD_unspecific|BD_I|BD_II|Schizoaffective_bipolar) and secondary diagnosis not allowed except for metabolic and eating disorders. Tauopathy is any individual with CERAD = {1}, Braak = {4,5,6} and secondary diagnosis allowed. PD is any individual with PD diagnosis (PD|PD_uncertain_plus_encephalitic) and secondary diagnosis can be only AD. FTD is any individual with FTD diagnosis (FTD) and secondary diagnosis can be only AD. All disease contrasts used in the study can be found in Supplementary Table 1.

Isolation and FANS of nuclei from frozen brain specimens with hashing

All buffers were supplemented with RNase inhibitors (Takara). Then, 25 mg of frozen postmortem human brain tissue was homogenized in cold lysis buffer (0.32 M sucrose, 5 mM CaCl₂, 3 mM magnesium acetate, 0.1 mM, EDTA, 10 mM Tris-HCl, pH8, 1 mM DTT, 0.1% Triton X-100) and filtered through a 40 µm cell strainer. The flow-through was underlaid with sucrose solution (1.8 M sucrose, 3 mM magnesium acetate, 1 mM DTT, 10 mM Tris-HCl, pH 8) and centrifuged at 107,000g for 1 h at 4 °C. The pellets were resuspended in PBS supplemented with 0.5% BSA. Six samples were processed in parallel. Up to 2 million nuclei from each sample were pelleted at 500g for 5 min at 4 °C. Nuclei were resuspended in 100 µl staining buffer (2% BSA, 0.02% Tween-20 in PBS) and incubated with 1 µg of a unique TotalSeq-A nuclear hashing antibody (BioLegend) for 30 min at 4 °C. Before fluorescence-activated nuclear sorting (FANS), volumes were brought up to 250 µl with PBS and 7-aminoactinomycin D (7-AAD) (Invitrogen) added according to the manufacturer’s instructions. 7-AAD-positive nuclei were sorted into tubes precoated with 5% BSA using the FACSAria flow cytometer (BD Biosciences) (Supplementary Fig. 30). FACSDiva software (BD Biosciences, v.8.0.2) was used for data collection.

snRNA-seq and hashing library preparation

After FANS, nuclei were subjected to two washes in 200 µl staining buffer, after which they were resuspended in 15 µl PBS and quantified (Countess II, Life Technologies). Concentrations were normalized and equal amounts of differentially hash-tagged nuclei were pooled. Using 10x Genomics single cell 3′ v3.1 reagents (10x Genomics), 60,000 (10,000 per donor) nuclei were run in each of two 10x Genomics lanes to create a technical replicate. At the cDNA amplification step (step 2.2) during library preparation, 1 µl of 2 µm HTO cDNA PCR ‘additive’ primer was added¹³². After cDNA amplification, supernatant from 0.6× SPRI selection was retained for HTO library generation. cDNA libraries were prepared according to the 10x Genomics protocol. HTO libraries were prepared as previously described¹³². cDNA and HTO libraries were sequenced at NYGC using the NovaSeq platform (Illumina).

Processing of snRNA-seq data

Alignment

Paired-end snRNA-seq library reads were aligned to the hg38 reference genome using STAR solo^133,134 and sample pools were demultiplexed using genotype matching via vireoSNP¹³⁵. After per-library count matrices were generated, cell by gene counts were aggregated into a single count matrix, which has 8,898,978 cells by 60,605 genes. The downstream processing was performed using Pegasus (v.1.7.0)¹³⁶ and scanpy (v.1.9.1)¹³⁷.

QC

We applied a rigorous three-step quality control (QC) process to filter out low-quality nuclei for subsequent downstream analyses. First, QC was applied at the cell level. Poor-quality nuclei were detected by thresholding based on UMI counts, gene counts and mitochondrial content. The QC thresholds were defined from log-transformed gene and UMI counts. Three median absolute deviations were used as the lower limit. This resulted in 986 for genes and 1,179 for UMI counts as lower bounds. We used a hard cut-off of 15,000 genes and 200,000 UMI counts as upper limits. Any cells with mitochondrial genes greater than 1% were filtered out. We checked for possible contamination from ambient RNA using CellBender¹³⁸. We also checked for cells with high fractions of reads that mapped to non-mRNAs, such as rRNA, sRNA and pseudogenes, as well as known confounding features, such as the lncRNA MALAT1. Second, QC was applied at the feature-level by removing those that were not robustly expressed in at least 0.05% of nuclei. Lastly, QC was applied at the donor-level by removing any sample represented by <50 nuclei. Finally, we excluded donors with low genotype concordance and sex discrepancies. Detailed QC statistics (median gene and UMI counts, median percent mitochondrial genes and cell counts), separated by 2,924 single-cell libraries used, before and after the removal of low-quality cells are available in Supplementary Table 12.

Normalization

We applied standard pegasus log-normalization with default parameters after QC.

Doublet detection

After QC and log-normalization, we performed the first-pass Leiden clustering analysis on the whole dataset to prepare input for doublet inference step. On the basis of the first-pass clusters, further filtering was carried out by removing inferred doublets using the Scrublet method¹³⁹.

Batch correction

We assessed the correlation between all pairs of donor and technical variables using canonical correlation analysis (Supplementary Fig. 7a). We attempted to remove unwanted technical batch effects during downstream analysis. During highly variable gene selection, we use the tissue source as a batch variable to remove the effect of brain tissue source. After the PCA step, the effects of UMI counts, mitochondrial percentage and the cell cycle were regressed out. We further removed the effect of sequencing batch using the Harmony method¹⁴⁰.

Defining cellular taxonomy using iterative clustering

Hierarchical cellular taxonomy

Cellular taxonomy was defined using a divide-and-conquer strategy. From the full dataset containing over 6 million nuclei, 8 major cell classes were defined using the following steps: 6,000 highly variable genes (HVGs) were selected from mean and dispersions trends¹⁴¹ using the default parameters (min_mean=0.0125, max_mean=3, min_disp=0.5) and the brain source as a batch variable after manually excluding sex and mitochondrial chromosomes. We used the k-nearest-neighbour (k-NN) graph calculated on the basis of harmony-corrected PCA embedding space to cluster nuclei of the same cell type using Leiden¹⁴² clustering algorithms. We used UMAP¹⁴³ to visualize the resulting clusters. From the class-level clusters, we subsetted the data by each class. Recalculating HVGs among cells in the same class enabled us to refocus on a feature space that is more relevant for the same class of cells. A k-NN graph was then calculated on the basis of the harmony-corrected PCA of the selected HVGs. Leiden clustering was used to annotate 27 subclass-level annotations. We iterated the same HVG–k-NN–Leiden clustering for all 27 subclasses yielding 67 subtypes of human brain cells. After obtaining annotations in three levels of hierarchy, the resulting clusters were aggregated into pseudobulk and cluster-wise Pearson correlation coefficient was calculated using existing human DLPFC¹⁰ and M1¹¹ annotations. We matched the annotations on the basis of both high correlation and specificity of the cell type. The final cellular taxonomy was compared to the following snRNA-seq datasets: PFC from ROSMAP cohort¹ and PFC (Brodmann area 9) from SEA-AD cohort² (Supplementary Table 11).

Evaluation of cellular taxonomy

To ensure that subsequent clustering truly identifies biological subtypes that are not confounded by technical variables such as brain sources or donor variables such as age and sex, we evaluated the quality of our most granular subtype-level cell annotations with other published datasets^1,2. To quantify the separability of cell groupings across datasets, we followed the approach of previous work in which local inverse Simpson’s index (LISI) was used to quantify the contribution of biological and technical factors to transcriptional variation¹⁴⁰. We computed the LISI scores using various covariates, including cell type, donor, age, tissue source and sex. LISI measures the degree of local label mixing in low-dimensional space and is sensitive to the number of distinct clusters present in a dataset. To ensure comparability across datasets, we applied LISI using a harmonized set of metadata categories and calculated scores using the compute_lisi function with a perplexity of 30 and nn_eps = 0, based on embeddings derived from a randomly sampled subset of cells. As datasets with more granular annotations (that is, greater numbers of clusters) tend to yield higher LISI values due to increased opportunities for label separation in local neighbourhoods, we normalized the mean LISI by dividing it by the natural logarithm of the number of clusters (adjusted LISI = mean LISI/log[n_clusters]). This log-based normalization compresses the scale of the cluster count, helping to avoid over-penalizing datasets with finer annotations while still accounting for granularity-related bias. The resulting adjusted LISI scores allow for more equitable comparisons of separability across datasets with varying label resolution.

Spatial validation of cellular taxonomy

Xenium in situ panel design

The Xenium Human Brain Gene Expression Panel (1000599, 10x Genomics) and a custom panel of 100 genes (Supplementary Table 10) were selected for the Xenium experiment. The 100 gene custom panel consisted mainly of subclass markers selected on the basis of specificity and gene expression level. The custom gene list was sent to 10x genomics, and the probe design was performed using their in-house pipeline.

Tissue preparation

Fresh-frozen tissue specimens of DLPFC were dissected into small blocks on ice. Tissue blocks were snap frozen by submerging in an isopentane (320404-1L, Sigma-Aldrich) bath chilled with dry ice and stored in −80 °C. Before cryosectioning, tissue blocks were allowed to equilibrate to the cryostat (HM505, Microm) chamber temperature, and were mounted with OCT (Tissue-Tek O.C.T. Compound, 4583, Sakura Finetek). After trimming, good-quality 10 µm sections were flattened on the cryostat stage and placed onto pre-equilibrated Xenium slides (Xenium Slides & Sample Prep Reagents, 1000460, 10x Genomics). Then, 2–3 sections were placed on each slide. The sections were further adhered by placing a finger on the backside of the slide for a few seconds and were then refrozen in the cryostat chamber. The slides were sealed in 50 ml tubes and stored at −80 °C until Xenium sample preparation.

Sample preparation

Xenium sample preparation was performed according to the manufacturer’s protocol (‘Xenium In Situ for Fresh Frozen Tissues – Fixation & Permeabilization, CG000581, Rev C’ and ‘Xenium In Situ Gene Expression - Probe Hybridization, Ligation & Amplification, User Guide, CG000582, Rev C’). In brief, fresh frozen sections mounted onto Xenium slides from the previous step were removed from −80 °C storage on dry ice before incubation at 37 °C for 1 min. The samples were then fixed in 4% paraformaldehyde (formaldehyde 16% in aqueous solution, 100503-917, VWR) in PBS for 30 min. After rinsing in PBS, the samples were permeabilized in 1% SDS for 2 min and then rinsed in PBS before being immersed in the prechilled 70% methanol and incubated for 60 min on ice. After rinsing the samples in PBS, the Xenium Cassettes were assembled onto the slides. The samples were incubated with a probe hybridization mix containing both the Xenium Human Brain Gene Expression Panel (1000599, 10x Genomics) and a 100 custom gene panel at 50 °C overnight to allow the probes to hybridize to targeted mRNAs. After probe hybridization, the samples were rinsed with PBST and incubated with Xenium post hybridization wash buffer at 37 °C for 30 min. The samples were then rinsed with PBST and ligation mix was added. Ligation was performed at 37 °C for 2 h to circularize the hybridized probes. After rinsing the samples with PBST, amplification master mix was added to enzymatically amplify the circularized probes at 30 °C for 2 h. After washing with TE buffer, autofluorescence was quenched according to the manufacturer’s protocol and nuclei were stained with DAPI before Xenium in situ analysis.

Data processing

The prepared samples were loaded into the Xenium analyser and run according to the manufacturer’s instructions ‘Xenium Analyzer User Guide CG000584 Rev B’. After the Xenium analyser was initiated, the correct gene panel was chosen, and decoding consumables (Xenium Decoding Consumables, PN-1000487, 10x Genomics) and reagents (Xenium Decoding Reagents, PN-1000461, 10x Genomics) were loaded. The bottom of the slides was carefully cleaned with ethanol before loading. Once the samples were loaded and the run was initiated, the instrument scanned the whole sample area of the slides using the DAPI channel, and regions of interest were selected to maximize the capture area. Results were generated by the instrument using the default settings. Instead of using 15 µm nuclei expansion distance for segmentation of cells, the default for the Xenium analyser, we applied nuclei only segmentation (Supplementary Fig. 4a) by resegmenting the results with 0 µm nuclei expansion using the Xenium ranger: xeniumranger resegment --id=demo --xenium-bundle=/path/to/xenium/files --expansion-distance=0 --resegment-nuclei=True.

QC and major cell type identification

After nuclei were segmented, cell × gene matrices were generated from the overlap of each segmented nuclear boundary with detected transcripts in the Xenium experiment. Nuclei were subsequently filtered by the number of detected transcripts (Supplementary Fig. 4b), and only those containing at least 30 nuclear transcripts were retained for downstream analysis (Supplementary Fig. 4c). Gene expression data from each sample were then log-normalized and normalized data were used for PCA, k-NN graph calculation and Leiden clustering. Clusters were then assigned to one of eight major cell types using label transfer (see the next section). Clusters in which <90% of nuclei shared the same predicted class, as well as nuclei whose predicted class differed from the cluster they were assigned to, were removed (Supplementary Fig. 4c). Finally, all retained nuclei were assigned the predicted class assigned to them (Supplementary Fig. 4d).

Class label prediction

We used scANVI¹⁴⁴ to perform reference-based label transfer from the RADC dataset. In brief, we performed the following steps. First, snRNA-seq gene expression data were subset to the genes shared with the Xenium gene panel. Next, we used the scvi-tools package^145,146 to train machine learning models for dimensionality reduction based on the reference dataset and its assigned labels (such as class and subclass). Models were run with 5 layers and 30 latent variables, and the scANVI model was trained for 20 epochs with a minimal sample of 100 cells per cluster per epoch. Lastly, a transfer model was trained for 100 epochs and applied to query data to assign labels on the basis of those the model was trained on from the reference. To assess the performance of each transfer model, we asked the model to predict labels in the reference data (using the subset gene pool) and evaluated the rate of correct prediction and biases in label misassignment for each predicted label.

Subclass and subtype label transfer for neuronal cells

After subsetting to nuclei labelled as EN or IN, these were then used as a query for a second scANVI label transfer—this one trained on subclass (or subtype) labels. Accuracy for both major cell type and subclass models assessed by predicted labels in the RADC dataset based only on the Xenium gene panel was estimated at >98%.

Pseudo-bulk concordance between subtypes in Xenium and RADC

To estimate the concordance of defined labels between the reference (RADC cohort) and Xenium datasets, we first generated pseudo-bulk estimates of each group (subtype) by summing the counts of each gene in all cells belonging to that group in each dataset. These pseudo-bulk count matrices were then subjected to TPM and log-transformation, before each gene was normalized by z scoring. The z-score matrices of the two datasets were then used to calculate the Pearson correlation between each pair of subtypes across the datasets.

Definition of rough laminar layers within tissue

After assignment of class and subclass identities to Xenium nuclei, we identified regional domains of interest and laminar layers using a thresholded local composition approach. First, we calculated the 30 nearest neighbours of each nucleus in each sample, using squidpy’s spatial_neighbors function¹⁴⁷ and constructed a matrix counting how many nearest neighbours of each subclass existed for each nucleus (also known as a neighbourhood composition vector¹⁴⁸). Next, we aggregated counts from EN subclasses corresponding to physically proximal laminar layers, to increase the signal. Specifically, we aggregated counts from EN_L3_5_IT_1, EN_L3_5_IT_2, EN_L3_5_IT_3 and EN_L5_ET into a single pseudolayer (L3–5), and similarly aggregated counts form EN_L6_IT_1, EN_L6_IT_2, EN_L6_IT_CT, EN_L6B and EN_L5_6_NP into another pseudolayer (L5–6). Counts from EN_L2_3_IT were treated as a laminar layer without aggregation (L2–3). After aggregation, we ran identical domain identification processes for four spatial domains, each with a corresponding nearest neighbour count column: white matter (oligodendrocytes), and each of the three laminar layers described above. Independent domain calling was performed as follows. (1) We generated an image representation of the relevant neighbour abundance count for each tissue section, with a pixel size of 25 μm. If multiple nuclei were assigned to the same pixel, neighbour abundances were averaged. Moreover, pixels with no assigned nuclei were assigned the value of their nearest non-empty pixel to avoid edge effects. (2) Next, the nearest neighbour count image was smoothed using Gaussian blurring (with a sigma of 4 pixels, equivalent to 100 μm) and smoothed values were associated with nuclei on the basis of the pixels they belonged to. An initial domain call was then generated by applying an Otsu threshold to the smoothed nearest neighbour counts per nucleus. (3) To further smooth the domain boundary and make it spatially continuous, we coarsened the image resolution and removed holes. In brief, the binary domain calls from step 2 were again used to generate a coarser image (with twice the pixel size = 50 μm). Pixel values were summarized as binary measures—any nucleus-containing pixels were assigned 1 if at least one of those nuclei was associated with the domain in the mask generated by step 2, and 0 otherwise—then empty pixels were assigned values on the basis of their nearest non-empty neighbours. Lastly, we applied scipy’s binary_fill_holes function to make the domain spatially continuous¹⁴⁹. (4) Values from this image were assigned to each associated nucleus to generate the final single-domain assignment. We applied this process to associate nuclei with each of the spatial domains described above independently and found that around 86.9% of nuclei were uniquely assigned to a single domain (~5.2% were not assigned a domain and ~7.8% were assigned to 2 or more).

Density estimation and comparison for IN subclasses

To estimate the density of each IN subclass within each spatial domain, we used the nucleus association from the laminar layer identification described above. To calculate density, we also needed to calculate the physical area occupied by each domain in each tissue section. This was done using the following method. (1) We generated a composite image with the same pixel size as the final domain calls (50 μm), in which pixels had a categorical value corresponding to one of the relevant domains if they were uniquely associated with that domain (and an NA value otherwise). (2) We then counted the number of pixels uniquely associated with each domain in each tissue section and converted it to μm². For density analysis, only nuclei uniquely associated with a single spatial domain were considered. Moreover, we excluded five samples from the comparative density analysis due to a lack of representation of at least one of the spatial domains identified above. Densities were then simply defined as the number of nuclei from each subclass uniquely associated with each of the three rough laminar domains, divided by the total area occupied by that domain in each relevant tissue section. Statistical testing between each collection of (n = 6) densities in physically proximal domains was then performed using a paired-sample Wilcoxon rank-sum test, as implemented by scipy’s Wilcoxon function¹⁴⁹. Density estimation and analysis for subtypes were performed similarly.

Processing of genotypes

DNA extraction and genotyping was performed as described previously¹⁵⁰. In brief, genomic DNA was extracted from frozen brain tissue using the QIAamp DNA Mini Kit (Qiagen), according to the manufacturer’s instructions. The samples were genotyped using the Infinium Psych Chip Array (Illumina) at the Mount Sinai Sequencing Core. Pre-imputation processing consisted of running the QC script HRC-1000G-check-bim.pl from the McCarthy Lab Group (https://www.well.ox.ac.uk/~wrayner/tools/), using the Trans-Omics for Precision Medicine (TOPMed)¹⁵¹. Genotypes were then phased and imputed on the TOPMed Imputation Server (https://imputation.biodatacatalyst.nhlbi.nih.gov). The samples with a mismatch between one’s self-reported and genetically inferred sex, suspected sex chromosome aneuploidies, high relatedness as defined by the KING kinship coefficient¹⁵² (KING > 0.177) and outlier heterozygosity (±3 s.d. from the mean) were removed. Moreover, samples with a sample-level missingness of >0.05 were removed and calculated within a subset of high-quality variants (variant-level missingness ≤ 0.02).

For ancestry assignment, genotypes were first merged with GRCh38 v2a 1000 Genomes Project data (https://wellcomeopenresearch.org/articles/4-50)¹⁵³ using BCFtools (v.1.9)¹⁵⁴. PLINK (v.2.0)¹⁵⁵ was then used to calculate the merged genotypes’ principal components (PCs), after filtering (minor allele frequency (MAF) ≥ 0.01, Hardy–Weinberg equilibrium (HWE) P ≥ 1 × 10⁻¹⁰, variant-level missingness ≤ 0.01, regions with high linkage disequilibrium (LD) removed) and LD pruning (window size = 1,000 kb, step size = 10, r² = 0.2) steps. For the samples of EUR ancestry assigned using the QDA method, autosomal biallelic variants with an imputation R² > 0.8, HWE P  ≥  1 × 10⁻⁶ and variant-level missingness ≤ 0.02 were retained. Genotypes were then annotated with ancestry-specific MAF values from the National Center for Biotechnology Information’s Allele Frequency Aggregator (ALFA) (https://ftp.ncbi.nih.gov/snp/population_frequency/latest_release/). Only variants with an ancestry-specific ALFA MAF ≥ 0.01 were retained.

PRS calculation

Polygenic risk scores (PRS) were computed for the PsychAD cohort using summary statistics from AD GWAS⁹⁶. The PRS-CS-auto method¹⁵⁶ was used, which incorporates continuous shrinkage priors to adjust the effect sizes from these summary statistics. An LD reference panel from the developers of PRS-CS, based on data from the 1000 Genomes Project¹⁵³, was used (https://github.com/getian107/PRScs). The default settings for PRS-CS were applied, including parameters a = 1 and b = 0.5 for the γ–γ prior, 1000 Markov chain Monte Carlo iterations, 500 burn-in iterations and a thinning factor of 5. The global shrinkage parameter phi was determined using a fully Bayesian method. PLINK (v.2.0)¹⁵⁵ was used to calculate the individual-level PRS. For AD PRS, two separate scores were calculated, both including and excluding the APOE locus. We conducted a sensitivity analysis comparing the PRS generated with and without the APOE locus. We found negligible differences between the two PRS, suggesting that the exclusion of APOE did not significantly alter the overall genetic risk profile (Supplementary Fig. 13b). The PRS with the APOE locus was used for the primary analysis.

Genetic heritability analysis of polygenic risk

We established a standardized pipeline for multi-marker analysis of GenoMic annotation (MAGMA) followed by scDRS. MAGMA incorporates the association P values of genetic variants from the latest genome-wide association study (GWAS). We used the following GWAS summary stats in the scDRS/MAGMA pipeline: AD⁹⁶ (excluding APOE locus), MS⁹⁷, PD⁹⁸, epilepsy⁹⁹, migraines¹⁰⁰, stroke¹⁰¹, ALS¹⁰², SCZ¹⁰³, BD¹⁰⁴, MDD¹⁰⁵, ASD¹⁰⁶, ADHD¹⁰⁷, insomnia¹⁰⁸, education¹⁰⁹, IQ¹¹⁰, alcoholism¹¹¹, OCD¹¹², Tourette syndrome¹¹³, obesity¹¹⁴, T2D¹¹⁵, cholesterol¹¹⁶, rheumatoid arthritis¹¹⁷, IBD¹¹⁸ and UC¹¹⁹. We applied MAGMA using a standard window of 35 kb upstream and 10 kb downstream around the gene body. We executed scDRS using the top 1,000 gene weights, sorted by z score. The MAGMA and scDRS pipeline were run using the following parameters. MAGMA was run using -snp-loc g1000_eur.bim (SNP location file corresponding to the Phase 3 1000 Genome Project) and --gene-loc NCBI38.gene.loc (gene location file from NCBI build 38). Both files were obtained from https://ctg.cncr.nl/software/magma. For scDRS, the default setting was used.

Statistical power and effect-size considerations in single-cell analyses

To complement the primary analyses presented in the Article, we conducted additional power and effect-size assessments to guide interpretation of the results and provide context for the sensitivity of our study. These analyses include both theoretical calculations and empirical evaluations designed to quantify the ability to detect changes in cellular composition and gene expression, as well as to highlight variability in power across cell types. On the basis of theoretical calculations, the statistical power to identify changes in cell type composition increases with the magnitude of the difference and the total number of cells observed. With the sample size and cell count in this study, we have >80% power to identify changes in cell type composition of ~0.1% (Supplementary Fig. 5a). For differential expression analyses, power increases with the sample size, standardized effect size and correlation between measured and true gene expression. Given the sample size in this study and moderated R² values between measured and true gene expression, we have >80% power to detect effect sizes of 1 on Cohen’s d scale (Supplementary Fig. 5b). In practice, there are many factors affecting statistical power in gene expression studies. In bulk RNA-seq experiments, power is determined by biological signal and sample size in addition to total read count and magnitude of gene expression. In single-cell datasets, additional technical factors include the number of cells collected for each cell cluster and the number of reads per cell. These numbers vary widely across cell types, as does the number of individuals with sufficient cell and read counts. Thus, empirical statistical power varies across cell types and genes even when testing the case/control analysis for one disease. Moreover, the Winner’s Curse phenomenon means that effect size estimates that pass a 5% FDR cut-off in low-powered cell types will be substantially overestimated compared with effect sizes in well-powered cell types. In our analysis of AD versus controls, there is wide variation in the number of expressed genes and the number of individuals passing filters based on cell and read count. We observe that the number of DEGs increases with the number of individuals passing filters and the number of expressed genes (Supplementary Fig. 5c). Yet the mean of the absolute value of the effect size estimates decreases with numbers of individuals and expressed genes (Supplementary Fig. 5d). Indeed, there is a strong negative correlation (Spearman’s rho = 0.58) between the number of DEGs and the mean absolute effect size (Supplementary Fig. 5e). Examining the estimated effect size in more detail, we see estimated effect size centred at zero for all cell subclasses (Supplementary Fig. 5f), but the mean of the absolute estimated effect sizes varies widely for DEGs (Supplementary Fig. 5g).

Variance partition analysis of gene expression

After creating a pseudobulk by aggregating single-cell library (Channel) per subclass-level annotation (Assay), expression values by assays were stacked using the StackedAssay function of Dreamlet. The resulting stacked pseudobulk enabled us to perform variance partition across cell types (stackedAssay). We used the following regression formula: Gene expression ~ (1|stackedAssay) + (1|Channel) + (1|SubID) + (1|Source) + (1|Ethnicity) + dx_bit + scale(Age) + Sex + scale(PMI) + log(n_genes) + percent_mito + mito_genes + ribo_genes + mito_ribo, where dx_bit indicates binary disease status excluding metabolic and eating disorders. Technical covariates log(n_genes), percent_mito, mito_genes, ribo_genes and mito_ribo were removed from the plotting and subsequent analysis because they explained less than 1 × 10⁻⁴% of overall gene expression variation.

To better interpret the variance partitioning analysis, we examined the relationship between variance explained and baseline gene expression across subclasses. Specifically, we plotted the proportion of variance attributed to donor, cell type and residual components as a function of mean expression (log₂[CPM]) (Supplementary Fig. 5h). Genes with higher expression showed greater variance explained by donor and cell type, and lower residual variance, reflecting a technical effect whereby higher read counts yield more precise estimates and therefore bias the analysis toward highly expressed genes.

We performed sensitivity analyses to ensure that performing variance partition analysis across all subclasses does not obscure the cell-type-specific signals. We performed variance partition analyses within more homogeneous cell populations by stacking cell-type-specific assays into three superclasses, namely, EN, IN and non-neuron, similar to previous approach by BICCN exploring interindividual variation of human brain¹⁵⁷. This setup enables better resolution of potentially meaningful covariates within a cell class while still adjusting for sample- and technical-level confounders. Overall, we observed highly concordant variance partition results between variance partition performed separately in each superclass and variance partition performed with all subclasses. The dominant source of variance remained the same as the CellType followed by the BrainDonor (Supplementary Fig. 6j). Spearman’s correlation of estimated variance between two models were high across all covariates, including brain donor, age, sex, ancestry and diagnosis (Supplementary Fig. 6k).

MAPT locus haplotyping

From our harmonized genotype calls, we selected common variants in 17q21.31 locus (chromosome 17: 45307631–46836264), performed PCA analysis of genotypes using 10 PCs, and used k-means clustering with k = 3 to call three genotype clusters, H1H1, H1H2 and H2H2. We additionally confirmed the H1 haplotype using two published SNPs, rs17763050 and rs8070723, known to associate³⁴. Haplotypes were estimated using Beagle (v.5.4)¹⁵⁸ on the selected genotypes of the 17q21.31 region. The estimation of the initial haplotype frequency model converged after one burn-in iteration, and the estimate of the genotype phase converged after 23 phasing iterations. For testing association with PD diagnosis, we used logistic regression with age, sex, 10 genotype PCs, and H1H1 status: PD ~ Age + Sex + Source + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + H1H1.

Moreover, we tested the contribution of the H1 haplotype to AD among non-APOE4 carriers³² but did not find a significant association (P ≤ 0.302). For testing association with AD diagnosis, we first subsetted for individuals who are not carriers of the APOE4 allele and tested for AD association using logistic regression with the formula: AD ~ Age + Sex + Source + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + H1H1.

Compositional variation analysis using crumblr

We applied the crumblr method (https://diseaseneurogenomics.github.io/crumblr) for testing the variation of cell type composition¹⁵⁹, as reflected in Figs. 4a–c and 6a–c. In summary, Crumblr scales the cell count ratio (that is, fractions) data using centred log-ratio (CLR) transformation and applies linear models. As CLR-transformed data are still highly heteroskedastic, the precision of measurements varies widely. Crumblr uses a fast asymptotic normal approximation of CLR-transformed counts from a Dirichlet-multinomial distribution to model the sampling variance of the transformed counts. Crumblr enables incorporating the sampling variance as precision weights to linear (mixed) models to increase power and control the false-positive rate. Crumblr also uses a variance stabilizing transform based on the precision weights to improve the performance of PCA and clustering. Hypothesis testing was computed using the following formula: Cell composition ~ scale(Age) + Sex + (phenotype of interest).

By including these variables, we account for potential confounders and improve the accuracy and reliability of our hypothesis testing.

Replication of compositional variation using SEA-AD MERFISH data

We downloaded MERFISH from the SEA-AD project as described in the original study². To compare results to our current work, we first created a slightly modified cell type annotation. In brief, we retained the subclass annotation present in the original data but merged some subclasses if multiple annotations in the MERFISH data mapped to a single subclass in the current study (Supplementary Table 11). We then performed compositional analysis using crumblr in the MERFISH data using this partially merged subclass annotation, using the following formula: Cell composition ~ (1|Donor.ID) + Sex + Age.at.Death + PMI + Brain.pH + (trait).

We performed this analysis using three AD-associated traits found both in the current study and the SEA-AD data: (1) CERAD; (2) Braak; and (3) cognitive status. We then compared the estimated log-transformed fold changes associated with each trait and each pair of matching subclasses (Supplementary Table 11) between the two studies. To jointly test the concordance of the compositional changes for all three traits, we transformed all log-transformed fold change values with the rank-based inverse normal transform, then tested for concordance using a linear mixed model with the lmerTest R package, using the following formula: RankNorm(logFC1) ~ RankNorm(logFC2) + (1|trait).

Validation of compositional variation using RNAscope

Tissue selection

We selected five AD and five neurotypical donors, matched for age and gender, and their DLPFC tissue blocks were obtained from the Mount Sinai Neuropathology Brain Bank and Research CoRE.

Probe design

We focused on four subclasses that showed robust variation in AD (Fig. 6a). Cell-type markers were selected on the basis of gene expression level and specificity (Fig. 2d and Supplementary Fig. 3i). The following probes were used: CUX2 (RNAscope probe Hs-CUX2, 425581) for EN_L2_3_IT, SST (RNAscope probe Hs-SST-C2, 310591-C2) for IN_SST, BMP5 (RNAscope probe Hs-BMP5-C3, 472461-C3) for VLMC and MYOCD (RNAscope probe Hs-MYOCD-C4, 416211-C4) and CASQ2 (RNAscope probe Hs-CASQ2-C3, 447111-C3) for SMCs.

Cryosectioning

Before sectioning, frozen brain tissue blocks were allowed to equilibrate to the cryostat (CryoStar NX50, Thermo Fisher Scientific) chamber temperature. Tissue blocks were mounted to chuck with OCT (Tissue-Tek O.C.T. Compound, 4583, Sakura Finetek). After trimming, regions of interest (ROIs) with high quality and intact grey and white matter were selected and scored with razor blades. Sections were cut at 16 µm thickness and collected onto the slides by gently touching the sections with the slides and two sections were collected per donor. The sections were further adhered to the slides by drying at room temperature before sealed in slide mailers (Fisherbrand 5-Place Slide Mailer, Thermo Fisher Scientific) and stored at −80 °C until the RNAscope assay.

RNAscope assay and TrueBlack treatment pipeline

Sample preparation was performed according to the ACD RNAscope Multiplex Fluorescent v2 for Fresh Frozen tissues protocol (UM 323100/Rev B) with an extended C2 TSA dye incubation of 45 min (determined on the basis of pilot experiment testing). The sections were fixed in freshly made 10% formalin (formaldehyde solution, F8775-25mL, Millipore Sigma) in 1× PBS at 4 °C for 1 h, after which the sections were rinsed twice with 1X PBS. The sections were dehydrated by placing the slides into a series of diluted ethanol and finally 100% ethanol. Before the creation of hydrophobic barriers, slides were air dried and thick hydrophobic barriers were drawn around each section with the hydrophobic barrier pen (ImmEDG Hydrophobic Barrier Pen, NC9545623, Vector Laboratories). After the barriers were dried, hydrogen peroxide was applied to the sections and incubated for 10 min at room temperature; after the removal of the hydrogen peroxide, the sections were washed twice with distilled water. The sections were then treated with protease IV for 30 min at room temperature and washed with distilled water. To hybridize the gene-specific probes, the sections were incubated with diluted and mixed probes for 2 h at 40 °C, then the sections were washed twice with 1× wash buffer for 2 min at room temperature. The sections were either stored in 5× SSC overnight or immediately processed with AMP hybridization. During AMP hybridization, the sections were incubated sequentially with AMP1, AMP2 and AMP3 for 30 min each at 40 °C, with two washes using 1× wash buffer for 2 min between each incubation. After the AMP steps, fluorescence signal development was performed by incubating with HRP-C1 for 15 min at 40 °C. After two 2 min washes with 1× wash buffer, a diluted fluorophore was added and incubated for 30 min at 40 °C, followed by two 2 min washes with 1× wash buffer. After the washes, HRP blocker was applied for 15 min at 40 °C, followed by two final 2 min washes with 1× wash buffer. The same procedure was repeated for HRP-C2 and HRP-C3, each with distinct fluorophores. Nuclei were counterstained with DAPI for 30 s, and washed with 1× wash buffer. After the RNAscope assay, the samples underwent a TrueBlack (20× TrueBlack reagent, 23007, Biotium) treatment according to the Biotium TrueBlack Lipofuscin Autofluorescence Quencher Protocol 2: Post-treatment with TrueBlack protocol (PSF006). Buffers on the sections were carefully removed, and 1× TrueBlack in 70% ethanol was applied to completely cover the tissue sections. After 30 s of incubation, the sections were washed three times with PBS. Coverslip mounting was performed using ProLong Gold Antifade Mountant (P36930, Thermo Fisher Scientific).

Imaging

RNAscope results were acquired with EVOS M7000 at ×40 (OLY XAPO 40X NA0.95/WD0.18, AMEP4907, Thermo Fisher Scientific). For each section, 2–3 ROIs were captured, the selected ROIs covered the entire grey matter and approximately two-thirds of the white matter, with good tissue quality and well represented the whole section. In total, four channels were imaged for nuclei and the C1–3 of the RNAscope signal, with DAPI, GFP, RFP and CY5 fluorophore cubes. At ×40, the ROIs contain 100–200 single fields of view, which were captured with the scan area function of the EVOS M7000. The scan was performed with serpentine horizontal scan pattern and quick scan model with less overlapping areas. Auto focus on every field of the DAPI channel was chosen to gain proper focus. Single-field images were saved as 16-bit TIF files.

Image processing and result analysis

Individual images were stitched into the whole ROIs and four channel images of the same ROI were stacked together using a custom Fiji ImageJ macro using the Grid/Collection Stitching Plugin. For the DAPI channel, the Grid: snake by rows function was used, and the compute_overlap and subpixel_accuracy options were chosen. After the DAPI channel images were stitched, the resulting ‘TileConfiguration.registered.txt’ files were used as guidance to stitch images of the other channels to make the resulting ROI images having the same dimension. Whole-ROI images of four channels were stacked using the Images to Stack and Stack to Hyperstack functions of Fiji imageJ. To calculate the gene expression level at the single-cell level, all stacked whole-ROI images were further analysed in quPATH, by counting the number of RNAscope punctate signals inside segmented cell borders. quPATH annotation was drawn to cover the whole ROI while avoiding areas of poor tissue quality or high-fluorescence background. The cells in the ROIs were segmented on the basis of DAPI nucleus signals with the cell detection function of the quPATH. Next, RNAscope punctate signals were detected by the subcellular detection function of the quPATH. The results were saved as csv files with the saveDetectionMeasurements function of the quPATH. To compare the difference in cell composition between the AD and control donors, a custom Python script was developed to quantify the cells of a certain cell type. The threshold of a positive cell is determined by evaluating the distribution of the cell and also by inspecting the images. The cell type was determined on the basis of the number of clusters detected for each channel (2 for CUX2, 0 for SST and BMP5). We observed a higher sensitivity for CUX2, so we used a higher minimum threshold to compensate. For each sample, cell counts and fraction were calculated. Cell fractions were transformed using the CLR method, which divides each part of the composition by the geometric mean of all parts and then takes the logarithm of the resulting ratio. Statistical differences between AD cases and controls (n = 10) were assessed using the Mann–Whitney–Wilcoxon test.

Differential gene expression analysis using Dreamlet

Owing to the increased variable complexity in a large-scale disease atlas, scaling single-cell based approaches to millions of cells across a wide range of phenotypes presents computational challenges¹⁶⁰ and can be suboptimal^{161,162,163,164}. To account for the scale of these data, complex study designs with repeated measures, and low read count per cell, we applied Dreamlet for differential expression analysis, which applies a pseudobulk approach, as reflected in Figs. 5a–c and 6d,e. Building from the previously developed statistical tool Dream¹⁶⁵, it applies linear mixed models to the differential expression problem in single-cell omics data. It starts by aggregating cells by the donor using a pseudobulk approach^161,162 and fits a regression model. For each feature and cell cluster, the following mixed model was applied: Gene expression ~ scale(Age) + Sex + scale(PMI) + log(n_genes) + percent_mito + mito_genes + mito_ribo + ribo_genes + (phenotype of interest), where categorical and numerical variables were modelled as random and fixed effects, respectively. If the phenotype of interest was a categorical variable, we set the intercept as 0 and used predefined contrasts between two factors. We performed a gene set analysis using the full spectrum of gene-level t-statistics using Zenith¹⁶⁶.

Note on the APOE genotype

We considered including the APOE genotype as a covariate in the differential expression analysis of AD, but this was not supported by the data. First, performing differential expression analysis on the basis of the APOE genotype did not identify any genes as differentially expressed at a study-wide FDR of 5%. Second, we then included the APOE genotype as a covariate in the differential expression analysis of AD versus controls. The estimated effect sizes for AD were very similar when comparing models with and without the APOE genotype (Supplementary Fig. 13a). No genes were significantly different in estimated effect size between the two models. On the basis of these empirical results, our analysis excluded the APOE genotype in the model.

Meta-analysis between brain sources

We conducted a meta-analysis to integrate results from different brain banks for the same disorder. Data tables from multiple brain banks were combined into a single list for each disorder and annotated with their respective sources. The meta_analysis() function in dreamlet was used to perform the meta-analysis, which involved combining data tables into a single data frame, grouping the data by assay. The s.e. was computed from the estimated log-transformed fold change and moderated t-statistics according to s.e. = log[FC]/t. The meta-analysis was performed using the rma() function from the metafor package¹⁶⁷ with a fixed-effects model. P values were adjusted using the FDR method, and the −log₁₀ of the FDR values were calculated. This method was applied to datasets of AD, DLBD, Vas, PD, Tau, FTD, SCZ and BD.

Meta-analysis across the same disease category

To further synthesize findings across multiple disorders, a meta-of-meta analysis was conducted, grouping the disorders into neurodegenerative and neuropsychiatric categories. The results from the initial meta-analyses for each disorder were combined into lists on the basis of their categories. The meta-of-meta function was used to perform this higher-level analysis. This function combined the meta-analysis results into a single data frame, grouped the data by assay and calculated the s.e. values using the formula abs(estimate/statistic). This approach was applied to create meta-of-meta analyses for all disorders, neurodegenerative disorders and neuropsychiatric disorders.

Evaluation of shared disease signatures

We modelled the total disease signature for a given condition as the sum of shared and disease-specific DEGs:

where Δτ_k denotes the total disease signature for a disease k, G_shared denotes a set of genes having cross-disorder effects, and G_distinct denotes a set of genes with non-shared (disease-specific) effects. This notation highlights that DEG profiles for each disease arise from contributions of both shared genes (common to all diseases) and distinct genes (specific to each disease).

To determine how DEG effects are shared across cell types and disease, we applied a multivariate Bayesian meta-analysis approach using the mashr software⁴⁸. The software uses a Bayesian approach to shrink effect sizes across genes and cell types to estimate the posterior effect sizes and posterior probability that an effect has the correct sign. In accompanying work, we have extended this approach to develop a formal statistical test to identify cell-type- and disease-specific effects¹⁶⁸. For this analysis, we conducted composite tests in a cell-type-specific manner, assessing the probability that a gene exhibits a non-zero effect across all eight disease contrasts. Genes with a posterior probability of ≥0.05 were classified as part of the shared component.

Construction of the correlation matrix

To calculate the correlation matrix, we used a systematic approach to quantify the relationships between genetic estimates across different neuropsychiatric and neurodegenerative disorders. Spearman’s correlation coefficients were calculated to assess the strength and direction of association between genetic estimates across different disorders. The calculation was performed for the common to each pair of disorders, grouped by assay, after the exclusion of the shared disease signatures.

The wide-format correlation matrix was converted to a matrix suitable for heat-map visualization. Missing values were replaced with zeros. Annotations indicating the number of significant genes were added to the rows and columns using the rowAnnotation and HeatmapAnnotation functions of the ComplexHeatmap package. The heat map was generated with hierarchical clustering of both rows and columns, and coloured on the basis of the Spearman’s correlation values using a gradient from blue (negative correlation) to red (positive correlation).

Co-heritability analysis

We used cross-trait LD score regression using the LDSC tool¹⁶⁹ to estimate the genetic correlation between a pair of traits. We used summary stats for the following GWAS traits (AD⁹⁶, DLBD¹⁷⁰, PD⁹⁸, FTD¹⁷¹, SCZ¹⁰³, BD¹⁰⁴) and calculated the heritability for each of the traits and the genetic covariance and correlation between each of the pair of traits (in total 15 pairs of traits). The size of the cohort was provided to the function munge_sumstats.py for heritability estimates. Precomputed LD scores for 1000 Genomes EUR data were downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/eur_w_ld_chr.tar.bz2. The SNP list for munge_sumstats.py was downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/w_hm3.snplist.bz2. Standard error was obtained from the LDSC output. The script munge_sumstats.py was modified to include the parameter --chunksize 5e5.

For each of the 15 possible combinations of traits, we calculated the level of correlation of gene expression using the Spearman rank correlation test. Genes were selected by applying the following criteria: log₂[FC] ≥ 0.5, FDR < 0.05. Co-expression coefficient was calculated for the overall dataset and for each of the cell types. Next, to correlate co-expression and co-heritability, we calculated the Spearman rank correlation coefficient between the LDSC genetic correlation score and the co-expression coefficient using 15 possible combinations of traits as datapoints for the Spearman rank correlation. Spearman was calculated for the overall dataset providing the genetic estimates of the expression similarities in the PsychAD cohort, and also per cell type to obtain the ranking of the cell types that contribute the most to the genetic to transcriptomic similarity in PsychAD.

Causal mediation analysis

Methods

Causal mediation analysis was performed on a subset of 645 individuals with European ancestry in AD contrast (tier 2, n = 696), who have PRS calculations from the latest AD GWAS⁹⁶. Two different R packages were used, namely mediation (https://cran.r-project.org/web/packages/mediation/) and psych (https://cran.r-project.org/web/packages/psych/). The results were cross-checked between the two methods (identical within a threshold) to ensure the estimated coefficients, and the mediation effects are statistically robust.

Statistical assumptions

A key assumption in causal mediation analysis is that there are no unobserved confounding effects between the mediator and the outcome. To address this, we carefully selected and included covariates in our model that are known or suspected to confound the mediator–outcome relationship. Furthermore, we ensured that the dependent, independent, and mediator variables were continuous, normally distributed and linearly related. For each regression, we used the following covariates: Age + Sex + PMI + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10, where PC1–10 indicate the genotype PCs. For the variation of subclass-level cell type composition, we used CLR-transformed cell count fractions from the crumblr compositional analysis. For bootstrapping, we used 10,000 simulations with the 50th percentiles of the treatment variable used as the control condition and the 90th percentile of the treatment variable used as the treatment condition, similar to a previous publication¹⁷².

Sensitivity analysis

We performed additional sensitivity analyses to strengthen the robustness of the mediation analysis. To ensure that the mediation effect was specific to AD, we performed the causal mediation analysis using PRSs derived from other neurological diseases (PD, MS, ALS, SCZ, BD and ASD). We found that none of these alternative PRSs yielded significant mediation effects (Supplementary Table 13), supporting the specificity of our findings to AD.

Trajectory analysis using neural network models

Rationale

Traditionally, changes in gene expression as a function of disease state are measured using linear-based models. This approach has proven to be highly valuable and has enhanced our understanding of the biological mechanisms underlying many diseases. However, it is increasingly recognized that changes in gene expression can be highly nonlinear⁶⁶; the interaction among numerous signalling pathways, many involving multiple feedback loops, can lead to complex dynamics that linear models may fail to capture. One approach to capturing potentially nonlinear changes in gene expression is pseudotime analysis (that is, trajectory inference)^173,174. Most trajectory inference tools rely on identifying continuous transitions within the dominant sources of variation in the data. However, as described in Fig. 3a, the variance in gene expression attributable to disease state is dwarfed by variance introduced by inter-donor heterogeneity, sex, brain source, and other technical or biological confounders. When disease related variance is not the principal axis of variation, these tools may construct trajectories that primarily reflect confounding factors rather than true disease progression. To address this limitation, we trained neural network models in a supervised manner that directly estimates the state of disease progression from single-cell gene expression data. In contrast to unsupervised trajectory inference methods, our approach explicitly models the relationship between gene expression and disease labels (for example, Braak stage and dementia status), enabling us to assign each cell a continuous disease progression score. This enables more accurate isolation and sensitive detection of disease-related signals, even when they account for a relatively small fraction of the overall variance.

Furthermore, in the case of AD, it is known that the spread of NFTs and Aβ plaques is strongly, but not perfectly, correlated with dementia. Of special interest is understanding cases in which individuals are resilient to dementia despite a high NFT or Aβ plaque burden. Thus, we aimed to disentangle changes in gene expression associated with disease burden from those associated with dementia. Standard trajectory inference methods do not allow us to separate these two covariates effectively. In our approach, we trained the neural network model by equally sampling all combinations of Braak stage and dementia status, thereby discouraging the model from learning spurious correlations between the two target variables. The model’s predictions of the Braak stage and dementia status were then used as two independent pseudotime axes.

Model architecture

The neural network model was a relatively simple feedforward network with two hidden layers. The input to the neural network were the log1p transformed gene counts from individual cells. Each hidden layer consisted of 1,024 units using the ReLU activation function. Layer and batch normalization were not used. The output of the network were Braak stage and the dementia status predictions.

Model training

Single-cell gene counts from the top 10,000 protein coding genes for each cell class, based on the percentage of cells the gene was expressed in, were used for training. We found that overfitting became more problematic when using a greater number of genes (data not shown). Genes found on sex chromosomes were excluded to discourage the model from learning to associate Braak or dementia status with sex. The network was trained to minimize the loss for both the Braak and dementia outputs. For the binary target dementia, the softmax function was applied to the output, and the loss was the cross-entropy. For Braak, target values were first normalized to zero mean and unit s.d., and the mean-squared error loss was used. Both loss terms were trained simultaneously. To prevent the model from overfitting the data, we applied dropout with probability of 0.5 to all hidden layers (after the ReLU activation). For each cell class, we divided the cells into 20 splits. In each split, ~95% of the cells was used for training, and the remaining ~5% was used for inference. Within each split, cells from a single donor exclusively belonged to either the training set or the inference set, but never both. Thus, model predictions were always based on cross-validated data from different donors. We trained one model for each of the 20 splits to generate predictions for all donors for that cell class. Models were trained using all donors from the MSSM and RADC brain banks and then ran inference on the 696 donors that focused on the AD phenotype contrast (Fig. 1d). There were different numbers of cells for each cell class; thus, the amount of training differed between classes. For neurons, astrocytes, oligodendrocytes and the immune cell class, we trained for 20 epochs. For mural and endothelial cells, which contained less data, we trained for 100 epochs. As model accuracy for the OPC cell class evolved more slowly over model training, we trained for 30 epochs. Cross-validated model accuracy across model training is shown in Supplementary Fig. 29. The models were trained using stochastic gradient descent (SGD) without momentum and with a learning rate of 0.02. We found that training with SGD led to greater accuracy compared with adaptive optimizers such as ADAM (data not shown). We used a batch size of 256 and linearly increased the learning rate across 2,000 training steps until it reached 0.02. Finally, we clipped the gradient norm to 1.0 to stabilize training. The model accuracy was determined as follows. We averaged the cell-level Braak and dementia model predictions to obtain donor-averaged scores. For dementia, we calculated the balanced classification score by determining the percentage of donors without dementia with a prediction score ≤ 0.5, the percentage of donors with dementia with a prediction score > 0.5 and then averaging these two values. For Braak, we calculated the Pearson R value between the actual Braak stage and the Braak predictions. Error bars were generated using a bootstrap procedure, in which we randomly sampled donors (with replacement), calculated the Braak and dementia prediction accuracy, and repeated this process 20,000 times. For each cell class, we included all donors with at least five cells. The number of donors with at least five cells and a defined Braak stage/Dementia status for each cell class was as follows: EN, 680/650; IN, 684/653; Astro, 692/662; Immune, 685/655; Oligo, 688/657; OPC, 688/659; Mural, 663/634; Endo, 589/560.

Disentangling Braak and dementia

Braak stage and dementia status are significantly correlated (Pearson R = 0.582, P < 1 × 10⁻⁶⁰). This strong correlation between target variables implies that input features (that is, changes in gene expression) associated with the two target variables are also likely to be correlated, which makes it challenging for the model to learn which input feature is predictive of which target variable. The result is when the model is trained in a standard manner, the Pearson correlation between the Braak and dementia model predictions is ~0.94 at the donor-level (Supplementary Fig. 14a; immune class shown), suggesting that the model has learned spurious correlations between the input features and target variables. Fully removing spurious correlation in machine learning models is still an unresolved question. However, balancing the training data, such that each of the 14 combinations of Braak and dementia (7 Braak values 2 dementia values) are equally sampled, can effectively reduce spurious correlations learned during training¹⁷⁵. In practice, we equally sampled from 15 groups, in which the extra group consisted of donors whose Braak stage or dementia status had not been determined. Training with group balancing reduced the correlation between the predicted Braak and dementia values (Supplementary Fig. 14b), did not adversely affect the model accuracy at the donor-level (Supplementary Fig. 14d).

Calculating gene trajectories

We wished to measure how gene expression varied as a function of the predicted Braak stage. First, gene counts were normalized so that each cell’s total count was 10,000, followed by the log1p transformation. Second, for each cell class and each donor, we calculated the mean predicted Braak stage (termed disease pseudotime), and the mean-normalized expression for each gene. Averaging within each donor reduced variability and ensured that donors with greater cell counts did not contribute disproportionately to downstream analysis. Third, we smoothed both the predicted donor-averaged disease pseudotime, and the donor-averaged gene expression with a Gaussian kernel. Specifically, for each donor i, we weighted all other donors j as

where p_i is the disease pseudotime of donor i, σ² was set to an eighth (1/8) of the variance of the disease pseudotime distribution and the normalization term Z_i was set such that \({\sum }_{j}{w}_{i,j}=1\). This enabled us to calculate smoothed disease pseudotime, \(\widetilde{{p}_{i}}\), and smoothed gene expression, \(\widetilde{{g}_{i}}\), values

where g_j is the gene expression vector of log1p normalized counts for donor j. After ordering the smoothed disease pseudotime values, gene trajectories are now represented as the tuple \((\widetilde{{p}_{i}},\widetilde{{g}_{i}})\). We only included donors with at least five cells for the cell class

Resilience against dementia

As tau proteinopathy and dementia status are highly correlated, gene expression as a function of the two variables is also correlated, and therefore partially redundant. We therefore aimed to measure how gene expression covaried with predicted dementia given the predicted Braak staging. To do so, we first calculated the expected predicted dementia and expected gene expression for donors with similar disease pseudotime. Specifically, we defined the expected dementia given disease pseudotime, \({d}_{i}^{p}\), and the expected gene expression given disease pseudotime, \({g}_{i}^{p}\),

where d_j is the predicted dementia status. As above, w_i,j was calculated on the basis of the difference in disease pseudotime between donors i and j, except that we set w_i,i = 0 so that each donor does not contribute to its own expected value. We then calculated the residuals between the donor’s predicted dementia status and gene expression with its expected values:

The dementia resilience score for each donor was then the product of these two terms. When calculating early and late resilience, donors were separated into early and late groups on the basis of disease pseudotime before averaging within each group. Using this metric, we define genes as protective if gene expression increases as predicted dementia decreases, given the disease pseudotime (that is, the product of the terms defined above is negative). Conversely, we define genes as damaging if gene expression increases as predicted dementia increases, given the disease pseudotime.

Identifying trajectory transition points

We sought to determine whether specific points during disease progression corresponded to shifts in gene dynamics. First, we included all donors with at least five cells in each of the eight cell classes (n = 578 donors). Next, we concatenated the smoothed trajectories of all coding genes across the eight classes (n = 17,265 genes × 8 cell classes) and performed PCA of the resulting 578 × 138,120 donor-by-gene/cell matrix. The trajectory projected onto the first three PCs, which together explain more than 91% of the variance, is shown in Supplementary Fig. 16, where the hue represents the disease pseudotime prediction averaged across all cell types for each donor. Visual inspection suggests points at which the trajectory changes direction. To quantify these transitions, we fit the trajectory projection (using the first three PCs) with a piecewise linear function consisting of two domains. Each domain represents a continuous block of donors ranked by mean disease pseudotime, separated by a transition point n (where n is the rank of the last donor in the first block). We calculated the explained variance for the piecewise fit across all possible transition points. The transition point that maximized explained variance was identified at n = 225 (out of 578 donors), which we define as the early-to-late transition point (Supplementary Fig. 16a (magenta circle)). A second, earlier transition is visible in the trajectory projection onto PC1 and PC3, although PC3 explains significantly less variance than PC2, making this shift weaker. To identify this preliminary transition, we repeated the same analysis using only the first 225 donors. This identified an earlier transition at n = 53 (out of 578 donors), marked by the cyan circle in Supplementary Fig. 16a. In subsequent analyses using these transition points to characterize dynamics within individual cell classes (Fig. 7), we adjusted the transition values proportionally on the basis of the number of donors available for each class. For example, in the EN cell class, 682 donors have at least five EN cells. To determine the early-to-late transition point for this class, we scaled the original transition point (225 out of 578 donors) according to the total number of EN donors: (682 × 225)/578 ≈ 265, rounding down.

Trajectory nonlinearity

We sought to determine whether gene trajectories were more nonlinear in certain cell classes than others. To assess this, we first fit the trajectory of all coding genes with a linear model and calculated the mean explained variance for each cell class (blue bars in Supplementary Fig. 17a). Next, we applied a piecewise linear fit, systematically varying the transition point across donors (as described in the ‘Identifying trajectory transition points’ section). For each cell class, we selected the transition point that maximized the mean explained variance across all coding genes. The explained variance at these optimal transition points is shown with orange bars in Supplementary Fig. 17a. To quantify trajectory nonlinearity, we defined a nonlinearity index as the difference between the explained variance of the optimal piecewise fit and that of the single linear fit, averaged across all coding genes (Supplementary Fig. 17b). Finally, we compared the nonlinearity index between neuronal and non-neuronal cell classes using a Wilcoxon rank-sum test.

Trajectory gene enrichment

We wished to determine which genetic pathways were most significantly up or downregulated during the progression of AD. To do so, we first extracted the slopes of the early and late linear fits for the Braak trajectories, and the mean early and late resilience scores (defined above). We used these slopes input to Zenith (https://bioconductor.org/packages/release/bioc/html/zenith.html) to calculate the changes across all GO BP pathways across the eight cell classes. For each GO BP pathway, we calculated the minimum Zenith FDR across the eight cell classes and across early and late stages, for both disease pseudotime and resilience. Only pathways with a minimum FDR ≤ 0.01 were included. Furthermore, we included only pathways with at least ten genes to ensure that the results were statistically robust, and no more than 250 genes to ensure that the pathways were not overly broad. Next, as we were interested only in pathways that could be informative of the mechanisms underlying AD progression, we excluded pathways containing words referring to overly broad behaviours or cognitive functions (learning, memory, vocalization, social, auditory, startle response, behaviour, locomotor, startle, prepulse inhibition), terms referring to anatomical structures other than the cortex (substantia nigra development, cardiac, coronary, aortic, ventricular, kidney, metanephric, retina, optic, bone, respiratory, pulmonary, olfactory, sperm, placenta, egg, embryonic, ovulation, estrous, placenta, sperm, mamary, germ layer, outflow tract septum, adrenal, epithelial, skeletal, otic, head) or overly broad neural terms (nervous system process, cerebral cortex, recognition, host, organ, developmental growth). To condense the remaining pathways into a more manageable size, we used rrvgo (https://www.bioconductor.org/packages/release/bioc/html/rrvgo.html). We selected the Wang semantic similarity metric¹⁷⁶, and set the threshold at 0.8 to obtain 86 GO BP pathways (Supplementary Fig. 18). For easier visualization, we selected 32 representative pathways from this set for Fig 7e. We also performed similar steps to obtain the top GO Molecular Function (MF) and Cellular Component (CC) pathways (Supplementary Figs. 19 and 20). As above, only pathways with a minimum FDR ≤ 0.01 were included, but the rrvgo threshold was set to 0.5 as there were fewer significant pathways.

Comparison to PLS regression

To compare our neural network approach with linear regression, we applied partial least squares (PLS) regression (using sklearn.cross_decomposition.PLSRegression) to predict Braak stage and dementia status in the immune cell class. As with our neural network models, the PLS model was jointly trained to predict Braak and dementia across 20 train/test splits, ensuring that predictions were generated on donors that were not included in model training. We empirically determined that setting n_components to 10 yielded the highest cross-validated accuracy. The PLS model achieved slightly greater Braak prediction accuracy than the neural network model, although the difference was not statistically significant (P > 0.05, bootstrap; Supplementary Fig. 15a). By contrast, the dementia classification accuracy was significantly lower (P < 0.001, bootstrap). However, we note that training a regression model on a binary classification task is not ideal, and improvements to this approach are possible. Moreover, PLS predictions for Braak and dementia were significantly more correlated than those from the neural network (Supplementary Fig. 15b). As Braak stage and dementia status are correlated, models that fail to disentangle their effects may overlook genes with opposing expression patterns. For example, if a gene’s expression increases with the Braak stage but decreases with dementia, these opposing trends could cancel out, rendering the change statistically insignificant. To evaluate this, we performed gene set enrichment analysis using Zenith on the early disease stage for both models, considering pathways for which the Zenith FDR was below 0.05 for either disease pseudotime or resilience. In Supplementary Fig. 15c, we compare the disease pseudotime Zenith z scores from the PLS and neural network models. The left panel shows pathways with concordant expression patterns (for example, increasing with Braak and associated with higher dementia risk), while the right panel highlights pathways with divergent patterns (for example, increasing with Braak but associated with lower dementia risk). For pathways with divergent trends, the PLS model yielded a lower Zenith z score magnitude. Finally, Supplementary Fig. 15d presents the top ten pathways with the greatest negative difference in z scores between the PLS and neural network models. These include key immune response pathways implicated in the early stages of AD.

Comparison to traditional trajectory inference methods

We wished to compare the accuracy of our neural network model predictions against two traditional trajectory inference methods: Monocle 3 (ref. ¹⁷⁷) and Palantir^177,178. For Monocle 3, we use the Python wrapper py-monocle (https://github.com/bioturing/py-monocle). We compared our method against these alternative approaches on the subset of immune cells consisting of microglia; we excluded adaptive and PVM cells as their different expression compared to microglia could be problematic for traditional trajectory inference approaches. The UMAP was computed as described in the ‘Defining cellular taxonomy using iterative clustering’ section. For both Monocle 3 and Palantir, we randomly selected a starting cell from among those with donor metadata indicating Braak stage = 0, CERAD score = 1, and no dementia or MCI. Trajectory inference was then performed, and the Pearson correlation between the resulting pseudotime and the actual Braak stage was calculated. This process was repeated 200 times, and the starting cell yielding the highest correlation was retained. The results in Supplementary Fig. 27a show the mean actual Braak, neural network model predicted Braak, Monocle3 pseudotime and Palantir pseudotime across UMAP space. The Pearson correlation between the actual Braak and model predictions or trajectory inference pseudotime values are shown in Supplementary Fig. 27b. Although it may be possible to optimize trajectory inference performance by refining hyperparameters, selecting different starting points or improving the underlying UMAP embedding, it is unclear whether such improvements would be sufficient to surpass the predictive accuracy of our neural network model (Supplementary Fig. 27b). Moreover, it is not evident how traditional trajectory inference approaches could be adapted to disentangle Braak stage from dementia status—an aspect our model addresses directly.

Pathway enrichment without smoothing

To ensure that our results were not biased by data smoothing (as described in the ‘Calculating gene trajectories’ section), we repeated our pathway enrichment analysis without smoothing gene trajectories. We maintained the same early-to-late transition point identified in Fig. 7c, which was derived from smoothed trajectories. To compute unsmoothed gene expression changes, we grouped donors into early and late stages on the basis of their mean disease pseudotime predictions. We then calculated slopes between donor-averaged gene expression and donor-averaged disease pseudotime predictions, without smoothing either variable. For resilience scores, we computed the partial covariance between donor-averaged gene expression and dementia predictions, regressing out disease pseudotime predictions. As before, these calculations were performed separately for early- and late-stage donors. Supplementary Fig. 28 compares Zenith z scores from the smoothed (x axis) and unsmoothed (y axis) approaches. Each dot represents a GO BP pathway, and we only included pathways in which the unsmoothed Zenith FDR was below 0.05. Rows correspond to cell classes, while columns compare early and late gene expression changes and resilience scores. Pearson correlations for increasing versus decreasing, or protective versus damaging pathways, ranged from 0.843 to 0.996 (median = 0.981). These results suggest that our findings are not an artifact of oversmoothing.

Identifying putative causal pathways

To identify GO BP pathways enriched in AD risk genes, we compared the z scores from the Alzheimer’s GWAS⁹⁶ between genes within each pathway and all other protein-coding genes not in the pathway. This was assessed using a one-sided t-test, and pathways with P values below 0.01 were considered to be potentially causal.

Mean trajectories and MAGMA enrichment

For both the mean normalized expression (Fig. 7f) and Magma enrichment analysis (Fig. 7h), we used the top 250 coding genes based on the early and late slopes of the Braak trajectories. Late-decreasing genes tend to also appear to contain an early increase (Fig. 7f). However, we cannot say whether this is biologically meaningful or the result of selection bias, as a strong late decrease must be preceded by a high baseline. For both the mean expression and MAGMA enrichment calculations, results were qualitatively similar if we used the top 500 or 1,000 genes instead (data not shown).

Immune cell dynamics

To further characterize how the immune response evolves during disease progression, we performed pathway enrichment using a sliding-window analysis. We sorted donors on the basis of their mean disease pseudotime predictions and calculated the change in pathway expression using the smoothed gene trajectories using windows of 60 donors and shifting the window by 20 donors at each interval. The Zenith z scores were calculated for each window of 60 donors. We included only pathways for which the Zenith FDRs were below 0.01 for at least 5 windows and included all pathways with between 8 and 250 genes. Furthermore, we only included putatively causal pathways (see the ‘Identifying putative causal pathways’ section), in addition to an immune and metabolic pathway that were presented in Fig. 7e (cytoplasmic translation and detection of bacterium). Using the mean disease pseudotime value from the midpoint of each window, we linearly interpolated the z-scores to obtain 100 values evenly spaced across the disease pseudotime axis. The vertical lines indicating the preliminary and early-to-late transition points were obtained from the PCA analysis in Supplementary Fig 16.

Ethics oversight

All procedures and research protocols were approved by the respective ethical committees of our collaborator’s institutions. The Ethics committee/IRB of Mount Sinai gave ethical approval for this work. The ethics committee/IRB of James J. Peters Department of Veterans Affairs Medical Center gave ethical approval for this work. The ethics committee/IRB of Rush Alzheimer’s Disease Center gave ethical approval for this work. The ethics committee/IRB of National Institute of Mental Health Human Brain Collection Core gave ethical approval for this work.

Reporting summary

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

All the source codes used in this study are available at GitHub (https://github.com/DiseaseNeuroGenomics/PsychADxD).

1. 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 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

2. Gabitto, M. I. et al. Integrated multimodal cell atlas of Alzheimer’s disease. Nat. Neurosci. 27, 2366–2383 (2024).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

3. Gandal, M. J. et al. Shared molecular neuropathology across major psychiatric disorders parallels polygenic overlap. Science 359, 693–697 (2018).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

4. Wingo, T. S. et al. Shared mechanisms across the major psychiatric and neurodegenerative diseases. Nat. Commun. 13, 4314 (2022).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

5. Mirra, S. S. et al. The Consortium to Establish a Registry for Alzheimer’s Disease (CERAD). Part II. Standardization of the neuropathologic assessment of Alzheimer’s disease. Neurology 41, 479–486 (1991).

   Article  CAS  PubMed  Google Scholar 

6. Braak, H. & Braak, E. Neuropathological stageing of Alzheimer-related changes. Acta Neuropathol. 82, 239–259 (1991).

   Article  CAS  PubMed  Google Scholar 

7. Braak, H., Alafuzoff, I., Arzberger, T., Kretzschmar, H. & Del Tredici, K. Staging of Alzheimer disease-associated neurofibrillary pathology using paraffin sections and immunocytochemistry. Acta Neuropathol. 112, 389–404 (2006).

   Article  PubMed  PubMed Central  Google Scholar 

8. Braak, H., Thal, D. R., Ghebremedhin, E. & Del Tredici, K. Stages of the pathologic process in Alzheimer disease: age categories from 1 to 100 years. J. Neuropathol. Exp. Neurol. 70, 960–969 (2011).

   Article  CAS  PubMed  Google Scholar 

9. Argyriou, S. et al. Beyond memory impairment: the complex phenotypic landscape of Alzheimer’s disease. Trends Mol. Med. https://doi.org/10.1016/j.molmed.2024.04.016 (2024).

10. Ma, S. et al. Molecular and cellular evolution of the primate dorsolateral prefrontal cortex. Science 377, eabo7257 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

11. BRAIN Initiative Cell Census Network (BICCN). A multimodal cell census and atlas of the mammalian primary motor cortex. Nature 598, 86–102 (2021).

    Article  Google Scholar 

12. Kwan, K. Y., Sestan, N. & Anton, E. S. Transcriptional co-regulation of neuronal migration and laminar identity in the neocortex. Development 139, 1535–1546 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

13. Markram, H. et al. Interneurons of the neocortical inhibitory system. Nat. Rev. Neurosci. 5, 793–807 (2004).

    Article  CAS  PubMed  Google Scholar 

14. Druga, R., Salaj, M. & Al-Redouan, A. Parvalbumin—positive neurons in the neocortex: a review. Physiol. Res. 72, S173–S191 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

15. Jagadeesh, K. A. et al. Identifying disease-critical cell types and cellular processes by integrating single-cell RNA-sequencing and human genetics. Nat. Genet. 54, 1479–1492 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

16. Zhang, M. J. et al. Polygenic enrichment distinguishes disease associations of individual cells in single-cell RNA-seq data. Nat. Genet. 54, 1572–1580 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

17. das Neves, R. P. et al. Connecting variability in global transcription rate to mitochondrial variability. PLoS Biol. 8, e1000560 (2010).

    Article  PubMed  PubMed Central  Google Scholar 

18. Johnston, I. G. et al. Mitochondrial variability as a source of extrinsic cellular noise. PLoS Comput. Biol. 8, e1002416 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

19. Zhu, X., Bührer, C. & Wellmann, S. Cold-inducible proteins CIRP and RBM3, a unique couple with activities far beyond the cold. Cell. Mol. Life Sci. 73, 3839–3859 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

20. Zhu, X. et al. Cold-inducible RNA binding protein alleviates iron overload-induced neural ferroptosis under perinatal hypoxia insult. Cell Death Differ. 31, 524–539 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

21. Ganat, Y., Soni, S., Chacon, M., Schwartz, M. L. & Vaccarino, F. M. Chronic hypoxia up-regulates fibroblast growth factor ligands in the perinatal brain and induces fibroblast growth factor-responsive radial glial cells in the sub-ependymal zone. Neuroscience 112, 977–991 (2002).

    Article  CAS  PubMed  Google Scholar 

22. Boettger, L. M., Handsaker, R. E., Zody, M. C. & McCarroll, S. A. Structural haplotypes and recent evolution of the human 17q21.31 region. Nat. Genet. 44, 881–885 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

23. Wade-Martins, R. Genetics: the MAPT locus-a genetic paradigm in disease susceptibility. Nat. Rev. Neurol. 8, 477–478 (2012).

    Article  CAS  PubMed  Google Scholar 

24. Tian, Y. et al. Shared genetics and comorbid genes of amyotrophic lateral sclerosis and Parkinson’s disease. Mov. Disord. 38, 1813–1821 (2023).

    Article  CAS  PubMed  Google Scholar 

25. Català-Senent, J. F. et al. A deep transcriptome meta-analysis reveals sex differences in multiple sclerosis. Neurobiol. Dis. 181, 106113 (2023).

    Article  PubMed  Google Scholar 

26. Cooper, Y. A. et al. Functional regulatory variants implicate distinct transcriptional networks in dementia. Science 377, eabi8654 (2022).

    Article  CAS  PubMed  Google Scholar 

27. Vandrovcova, J. et al. Association of MAPT haplotype-tagging SNPs with sporadic Parkinson’s disease. Neurobiol. Aging 30, 1477–1482 (2009).

    Article  CAS  PubMed  Google Scholar 

28. Steinberg, K. M. et al. Structural diversity and African origin of the 17q21.31 inversion polymorphism. Nat. Genet. 44, 872–880 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Stefansson, H. et al. A common inversion under selection in Europeans. Nat. Genet. 37, 129–137 (2005).

    Article  CAS  PubMed  Google Scholar 

30. Healy, D. G. et al. Tau gene and Parkinson’s disease: a case-control study and meta-analysis. J. Neurol. Neurosurg. Psychiatry 75, 962–965 (2004).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

31. Zabetian, C. P. et al. Association analysis of MAPT H1 haplotype and subhaplotypes in Parkinson’s disease. Ann. Neurol. 62, 137–144 (2007).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

32. Sánchez-Juan, P. et al. The MAPT H1 haplotype is a risk factor for Alzheimer’s disease in APOE ε4 non-carriers. Front. Aging Neurosci. 11, 327 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

33. Koks, S., Pfaff, A. L., Bubb, V. J. & Quinn, J. P. Transcript variants of genes involved in neurodegeneration are differentially regulated by the APOE and MAPT haplotypes. Genes 12, 423 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

34. Bowles, K. R. et al. 17q21.31 sub-haplotypes underlying H1-associated risk for Parkinson’s disease are associated with LRRC37A/2 expression in astrocytes. Mol. Neurodegener. 17, 48 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

35. Edwards, T. L. et al. Genome-wide association study confirms SNPs in SNCA and the MAPT region as common risk factors for Parkinson disease. Ann. Hum. Genet. 74, 97–109 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

36. Hill-Burns, E. M. et al. Identification of a novel Parkinson’s disease locus via stratified genome-wide association study. BMC Genom. 15, 118 (2014).

    Article  Google Scholar 

37. Erro, R., Mencacci, N. E. & Bhatia, K. P. The emerging role of phosphodiesterases in movement disorders. Mov. Disord. 36, 2225–2243 (2021).

    Article  PubMed  PubMed Central  Google Scholar 

38. Appenzeller, S. et al. Autosomal-dominant striatal degeneration is caused by a mutation in the phosphodiesterase 8B gene. Am. J. Hum. Genet. 86, 83–87 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

39. Azuma, R. et al. A novel mutation of PDE8B gene in a Japanese family with autosomal-dominant striatal degeneration. Mov. Disord. 30, 1964–1967 (2015).

    Article  CAS  PubMed  Google Scholar 

40. Barsottini, O. G. P. et al. Familial striatal degeneration: new mutation and neuroimaging clues. Neurology 85, 1816–1818 (2015).

    Article  PubMed  Google Scholar 

41. Folci, A. et al. Loss of hnRNP K impairs synaptic plasticity in hippocampal neurons. J. Neurosci. 34, 9088–9095 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

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

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

43. Zeng, T., Spence, J. P., Mostafavi, H. & Pritchard, J. K. Bayesian estimation of gene constraint from an evolutionary model with gene features. Nat. Genet. 56, 1632–1643 (2024).

    Article  CAS  PubMed  Google Scholar 

44. Jovanovic, M. et al. Dynamic profiling of the protein life cycle in response to pathogens. Science 347, 1259038 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

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

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

46. Yang, A. C. et al. A human brain vascular atlas reveals diverse mediators of Alzheimer’s risk. Nature 603, 885–892 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

47. 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).

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

    Article  CAS  PubMed  Google Scholar 

49. Boyle, E. A., Li, Y. I. & Pritchard, J. K. An expanded view of complex traits: from polygenic to omnigenic. Cell 169, 1177–1186 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

50. Weiner, D. J. et al. Polygenic architecture of rare coding variation across 394,783 exomes. Nature 614, 492–499 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

51. Sadeghi, I. et al. Brain transcriptomic profiling reveals common alterations across neurodegenerative and psychiatric disorders. Comput. Struct. Biotechnol. J. 20, 4549–4561 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

52. Voineagu, I. et al. Transcriptomic analysis of autistic brain reveals convergent molecular pathology. Nature 474, 380–384 (2011).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

53. Lyketsos, C. G. et al. Prevalence of neuropsychiatric symptoms in dementia and mild cognitive impairment: results from the cardiovascular health study. JAMA 288, 1475–1483 (2002).

    Article  PubMed  Google Scholar 

54. Bhattacherjee, A. et al. Cell type-specific transcriptional programs in mouse prefrontal cortex during adolescence and addiction. Nat. Commun. 10, 4169 (2019).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

55. Green, G. S. et al. Cellular communities reveal trajectories of brain ageing and Alzheimer’s disease. Nature 633, 634–645 (2024).

56. Kosoy, R. et al. Alzheimer's disease transcriptional landscape in ex vivo human microglia. Nat. Neurosci. 28, 1830–1843 (2025).

57. Lee, D. et al. Plasticity of human microglia and brain perivascular macrophages in aging and Alzheimer’s disease. Nat. Genet. https://doi.org/10.1038/s41588-026-02716-6 (2026).

58. Galea, I. The blood-brain barrier in systemic infection and inflammation. Cell. Mol. Immunol. 18, 2489–2501 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

59. Heneka, M. T. et al. Neuroinflammation in Alzheimer disease. Nat. Rev. Immunol. https://doi.org/10.1038/s41577-024-01104-7 (2024).

60. Fan, Z., Brooks, D. J., Okello, A. & Edison, P. An early and late peak in microglial activation in Alzheimer’s disease trajectory. Brain 140, 792–803 (2017).

    PubMed  PubMed Central  Google Scholar 

61. Haney, M. S. et al. APOE4/4 is linked to damaging lipid droplets in Alzheimer’s disease microglia. Nature 628, 154–161 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

62. Marschallinger, J. et al. Lipid-droplet-accumulating microglia represent a dysfunctional and proinflammatory state in the aging brain. Nat. Neurosci. 23, 194–208 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

63. Claes, C. et al. Plaque-associated human microglia accumulate lipid droplets in a chimeric model of Alzheimer’s disease. Mol. Neurodegener. 16, 50 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

64. Grayson, J. M. et al. T cell exhaustion is associated with cognitive status and amyloid accumulation in Alzheimer’s disease. Sci. Rep. 13, 15779 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

65. Chen, J., Liu, X. & Zhong, Y. Interleukin-17A: the key cytokine in neurodegenerative diseases. Front. Aging Neurosci. 12, 566922 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

66. Kitano, H. Computational systems biology. Nature 420, 206–210 (2002).

    Article  ADS  CAS  PubMed  Google Scholar 

67. Bhembre, N., Bonthron, C. & Opazo, P. Synaptic compensatory plasticity in Alzheimer’s disease. J. Neurosci. 43, 6833–6840 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

68. Hernández-Ortega, K., Garcia-Esparcia, P., Gil, L., Lucas, J. J. & Ferrer, I. Altered machinery of protein synthesis in Alzheimer’s: from the nucleolus to the ribosome. Brain Pathol. 26, 593–605 (2016).

    Article  PubMed  Google Scholar 

69. Moreira, P. I., Carvalho, C., Zhu, X., Smith, M. A. & Perry, G. Mitochondrial dysfunction is a trigger of Alzheimer’s disease pathophysiology. Biochim. Biophys. Acta 1802, 2–10 (2010).

    Article  CAS  PubMed  Google Scholar 

70. Lee, J.-H. et al. Faulty autolysosome acidification in Alzheimer’s disease mouse models induces autophagic build-up of Aβ in neurons, yielding senile plaques. Nat. Neurosci. 25, 688–701 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

71. Butterfield, D. A. & Halliwell, B. Oxidative stress, dysfunctional glucose metabolism and Alzheimer disease. Nat. Rev. Neurosci. 20, 148–160 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

72. Stefanoska, K. et al. Alzheimer’s disease: ablating single master site abolishes tau hyperphosphorylation. Sci. Adv. https://doi.org/10.1126/sciadv.abl8809 (2022).

73. Schwartzentruber, J. et al. Genome-wide meta-analysis, fine-mapping and integrative prioritization implicate new Alzheimer’s disease risk genes. Nat. Genet. 53, 392–402 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

74. Mallach, A. et al. Microglia-astrocyte crosstalk in the amyloid plaque niche of an Alzheimer’s disease mouse model, as revealed by spatial transcriptomics. Cell Rep. 43, 114216 (2024).

    Article  CAS  PubMed  Google Scholar 

75. Marsh, S. E. et al. The adaptive immune system restrains Alzheimer’s disease pathogenesis by modulating microglial function. Proc. Natl Acad. Sci. USA 113, E1316–E1325 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

76. Su, W. et al. CXCR6 orchestrates brain CD8⁺ T cell residency and limits mouse Alzheimer’s disease pathology. Nat. Immunol. 24, 1735–1747 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

77. Louveau, A. et al. Structural and functional features of central nervous system lymphatic vessels. Nature 523, 337–341 (2015).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

78. Ahn, J. H. et al. Meningeal lymphatic vessels at the skull base drain cerebrospinal fluid. Nature 572, 62–66 (2019).

    Article  ADS  CAS  PubMed  Google Scholar 

79. Li, X. et al. Meningeal lymphatic vessels mediate neurotropic viral drainage from the central nervous system. Nat. Neurosci. 25, 577–587 (2022).

    Article  CAS  PubMed  Google Scholar 

80. Fitzpatrick, Z. et al. Venous-plexus-associated lymphoid hubs support meningeal humoral immunity. Nature 628, 612–619 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

81. Da Mesquita, S. et al. Meningeal lymphatics affect microglia responses and anti-Aβ immunotherapy. Nature 593, 255–260 (2021).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

82. Louveau, A. et al. CNS lymphatic drainage and neuroinflammation are regulated by meningeal lymphatic vasculature. Nat. Neurosci. 21, 1380–1391 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

83. Da Mesquita, S., Fu, Z. & Kipnis, J. The meningeal lymphatic system: a new player in neurophysiology. Neuron 100, 375–388 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

84. Da Mesquita, S. et al. Functional aspects of meningeal lymphatics in ageing and Alzheimer’s disease. Nature 560, 185–191 (2018).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

85. Long, J. M. & Holtzman, D. M. Alzheimer disease: an update on pathobiology and treatment strategies. Cell 179, 312–339 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

86. Ryder, B. D., Wydorski, P. M., Hou, Z. & Joachimiak, L. A. Chaperoning shape-shifting tau in disease. Trends Biochem. Sci. 47, 301–313 (2022).

    Article  CAS  PubMed  Google Scholar 

87. Mok, S.-A. et al. Mapping interactions with the chaperone network reveals factors that protect against tau aggregation. Nat. Struct. Mol. Biol. 25, 384–393 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

88. Koren, J. III et al. Chaperone signalling complexes in Alzheimer’s disease. J. Cell. Mol. Med. 13, 619–630 (2009).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

89. Saha, I. et al. The AAA+ chaperone VCP disaggregates Tau fibrils and generates aggregate seeds in a cellular system. Nat. Commun. 14, 560 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

90. Abellanas, M. A., Purnapatre, M., Burgaletto, C. & Schwartz, M. Monocyte-derived macrophages act as reinforcements when microglia fall short in Alzheimer’s disease. Nat. Neurosci. 28, 436–445 (2025).

    Article  CAS  PubMed  Google Scholar 

91. Brigas, H. C. et al. IL-17 triggers the onset of cognitive and synaptic deficits in early stages of Alzheimer’s disease. Cell Rep. 36, 109574 (2021).

    Article  CAS  PubMed  Google Scholar 

92. Rosenzweig, N. et al. Sex-dependent APOE4 neutrophil-microglia interactions drive cognitive impairment in Alzheimer’s disease. Nat. Med. https://doi.org/10.1038/s41591-024-03122-3 (2024).

93. Cipollini, V., Anrather, J., Orzi, F. & Iadecola, C. Th17 and cognitive impairment: possible mechanisms of action. Front. Neuroanat. 13, 95 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

94. Criado-Marrero, M. et al. Hsp90 co-chaperones, FKBP52 and Aha1, promote tau pathogenesis in aged wild-type mice. Acta Neuropathol. Commun 9, 65 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

95. Gorantla, N. V. & Chinnathambi, S. Tau protein squired by molecular chaperones during Alzheimer’s disease. J. Mol. Neurosci. 66, 356–368 (2018).

    Article  CAS  PubMed  Google Scholar 

96. Bellenguez, C. et al. New insights into the genetic etiology of Alzheimer’s disease and related dementias. Nat. Genet. 54, 412–436 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

97. International Multiple Sclerosis Genetics Consortium. Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science 365, eaav7188 (2019).

    Article  Google Scholar 

98. 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).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

99. International League Against Epilepsy Consortium on Complex Epilepsies. Genetic determinants of common epilepsies: a meta-analysis of genome-wide association studies. Lancet Neurol. 13, 893–903 (2014).

    Article  Google Scholar 

100. 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).

     Article  PubMed  PubMed Central  Google Scholar 

101. Malik, R. et al. Genome-wide meta-analysis identifies 3 novel loci associated with stroke. Ann. Neurol. 84, 934–939 (2018).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

102. 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).

     Article  PubMed  PubMed Central  Google Scholar 

103. Trubetskoy, V. et al. Mapping genomic loci implicates genes and synaptic biology in schizophrenia. Nature 604, 502–508 (2022).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

104. 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).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

105. 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).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

106. Grove, J. et al. Identification of common genetic risk variants for autism spectrum disorder. Nat. Genet. 51, 431–444 (2019).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

107. 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).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

108. 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).

     Article  CAS  PubMed  Google Scholar 

109. Lee, J. J. et al. Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals. Nat. Genet. 50, 1112–1121 (2018).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

110. Sniekers, S. et al. Genome-wide association meta-analysis of 78,308 individuals identifies new loci and genes influencing human intelligence. Nat. Genet. 49, 1107–1112 (2017).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

111. 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).

     Article  PubMed  Google Scholar 

112. 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).

     Article  PubMed  Google Scholar 

113. Yu, D. et al. Interrogating the genetic determinants of Tourette’s syndrome and other tic disorders through genome-wide association studies. Am. J. Psychiatry 176, 217–227 (2019).

     Article  PubMed  PubMed Central  Google Scholar 

114. Watanabe, K. et al. A global overview of pleiotropy and genetic architecture in complex traits. Nat. Genet. 51, 1339–1348 (2019).

     Article  CAS  PubMed  Google Scholar 

115. Morris, A. P. et al. Large-scale association analysis provides insights into the genetic architecture and pathophysiology of type 2 diabetes. Nat. Genet. 44, 981–990 (2012).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

116. Willer, C. J. et al. Discovery and refinement of loci associated with lipid levels. Nat. Genet. 45, 1274–1283 (2013).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

117. Okada, Y. et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature 506, 376–381 (2014).

     Article  ADS  CAS  PubMed  Google Scholar 

118. Liu, J. Z. et al. Association analyses identify 38 susceptibility loci for inflammatory bowel disease and highlight shared genetic risk across populations. Nat. Genet. 47, 979–986 (2015).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

119. Anderson, C. A. et al. Meta-analysis identifies 29 additional ulcerative colitis risk loci, increasing the number of confirmed associations to 47. Nat. Genet. 43, 246–252 (2011).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

120. Hsiao, L. L. et al. A compendium of gene expression in normal human tissues. Physiol. Genom. 7, 97–104 (2001).

     Article  CAS  Google Scholar 

121. Hounkpe, B. W., Chenou, F., de Lima, F. & De Paula, E. V. HRT Atlas v1.0 database: redefining human and mouse housekeeping genes and candidate reference transcripts by mining massive RNA-seq datasets. Nucleic Acids Res. 49, D947–D955 (2021).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

122. The International League Against Epilepsy Consortium on Complex Epilepsies. Genome-wide mega-analysis identifies 16 loci and highlights diverse biological mechanisms in the common epilepsies. Nat. Commun. 9, 5269 (2018).

123. Savage, J. E. et al. Genome-wide association meta-analysis in 269,867 individuals identifies new genetic and functional links to intelligence. Nat. Genet. 50, 912–919 (2018).

124. Bennett, D. A. et al. Religious Orders Study and Rush Memory and Aging Project. J. Alzheimers Dis. 64, S161–S189 (2018).

     Article  PubMed  PubMed Central  Google Scholar 

125. Marquez, D. X. et al. Representation of older Latinxs in cohort studies at the Rush Alzheimer’s Disease Center. Neuroepidemiology 54, 404–418 (2020).

     Article  PubMed  PubMed Central  Google Scholar 

126. Bennett, D. A. et al. Neuropathology of older persons without cognitive impairment from two community-based studies. Neurology 66, 1837–1844 (2006).

     Article  CAS  PubMed  Google Scholar 

127. Bennett, D. A. et al. Natural history of mild cognitive impairment in older persons. Neurology 59, 198–205 (2002).

     Article  CAS  PubMed  Google Scholar 

128. Bennett, D. A. et al. Decision rules guiding the clinical diagnosis of Alzheimer’s disease in two community-based cohort studies compared to standard practice in a clinic-based cohort study. Neuroepidemiology 27, 169–176 (2006).

     Article  PubMed  Google Scholar 

129. Fullard, J. F. et al. Population-scale cross-disorder atlas of the human prefrontal cortex at single-cell resolution. Sci. Data 12, 954 (2025).

130. Pedregosa, F. et al. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830 (2011).

     Google Scholar 

131. Borrell, L. N. et al. Race and genetic ancestry in medicine—a time for reckoning with racism. N. Engl. J. Med. 384, 474–480 (2021).

     Article  PubMed  PubMed Central  Google Scholar 

132. Stoeckius, M. et al. Cell hashing with barcoded antibodies enables multiplexing and doublet detection for single cell genomics. Genome Biol. 19, 224 (2018).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

133. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

     Article  CAS  PubMed  Google Scholar 

134. Kaminow, B., Yunusov, D. & Dobin, A. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data. Preprint at bioRxiv https://doi.org/10.1101/2021.05.05.442755 (2021).

135. Huang, Y., McCarthy, D. J. & Stegle, O. Vireo: Bayesian demultiplexing of pooled single-cell RNA-seq data without genotype reference. Genome Biol. 20, 273 (2019).

     Article  PubMed  PubMed Central  Google Scholar 

136. Li, B. et al. Cumulus provides cloud-based data analysis for large-scale single-cell and single-nucleus RNA-seq. Nat. Methods 17, 793–798 (2020).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

137. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

     Article  PubMed  PubMed Central  Google Scholar 

138. Fleming, S. J. et al. Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nat. Methods 20, 1323–1335 (2023).

     Article  CAS  PubMed  Google Scholar 

139. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281–291 (2019).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

140. Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16, 1289–1296 (2019).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

141. Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

142. Traag, V. A., Waltman, L. & van Eck, N. J. From Louvain to Leiden: guaranteeing well-connected communities. Sci. Rep. 9, 5233 (2019).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

143. McInnes, L., Healy, J. & Melville, J. UMAP: uniform manifold approximation and projection for dimension reduction. Preprint at https://arxiv.org/abs/1802.03426 (2018).

144. Xu, C. et al. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models. Mol. Syst. Biol. 17, e9620 (2021).

     Article  PubMed  PubMed Central  Google Scholar 

145. Gayoso, A. et al. A Python library for probabilistic analysis of single-cell omics data. Nat. Biotechnol. 40, 163–166 (2022).

     Article  CAS  PubMed  Google Scholar 

146. Virshup, I. et al. The scverse project provides a computational ecosystem for single-cell omics data analysis. Nat. Biotechnol. 41, 604–606 (2023).

     Article  CAS  PubMed  Google Scholar 

147. Palla, G. et al. Squidpy: a scalable framework for spatial omics analysis. Nat. Methods 19, 171–178 (2022).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

148. Petukhov, V. et al. Cell segmentation in imaging-based spatial transcriptomics. Nat. Biotechnol. 40, 345–354 (2022).

     Article  CAS  PubMed  Google Scholar 

149. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272 (2020).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

150. Kosoy, R. et al. Genetics of the human microglia regulome refines Alzheimer’s disease risk loci. Nat. Genet. 54, 1145–1154 (2022).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

151. Taliun, D. et al. Sequencing of 53,831 diverse genomes from the NHLBI TOPMed Program. Nature 590, 290–299 (2021).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

152. Manichaikul, A. et al. Robust relationship inference in genome-wide association studies. Bioinformatics 26, 2867–2873 (2010).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

153. 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature 526, 68–74 (2015).

     Article  Google Scholar 

154. Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10, giab008 (2021).

     Article  PubMed  PubMed Central  Google Scholar 

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

     Article  PubMed  PubMed Central  Google Scholar 

156. Ge, T., Chen, C.-Y., Ni, Y., Feng, Y.-C. A. & Smoller, J. W. Polygenic prediction via Bayesian regression and continuous shrinkage priors. Nat. Commun. 10, 1776 (2019).

     Article  ADS  PubMed  PubMed Central  Google Scholar 

157. Johansen, N. et al. Interindividual variation in human cortical cell type abundance and expression. Science 382, eadf2359 (2023).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

158. Browning, B. L., Tian, X., Zhou, Y. & Browning, S. R. Fast two-stage phasing of large-scale sequence data. Am. J. Hum. Genet. 108, 1880–1890 (2021).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

159. Hoffman, G. E. GabrielHoffman/Crumblr: freeze to Zenodo. Zenodo https://doi.org/10.5281/zenodo.12752107 (2024).

160. Lähnemann, D. et al. Eleven grand challenges in single-cell data science. Genome Biol. 21, 31 (2020).

     Article  PubMed  PubMed Central  Google Scholar 

161. Crowell, H. L. et al. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat. Commun. 11, 6077 (2020).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

162. Squair, J. W. et al. Confronting false discoveries in single-cell differential expression. Nat. Commun. 12, 5692 (2021).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

163. Zimmerman, K. D., Espeland, M. A. & Langefeld, C. D. A practical solution to pseudoreplication bias in single-cell studies. Nat. Commun. 12, 738 (2021).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

164. Murphy, A. E., Fancy, N. & Skene, N. G. Avoiding false discoveries in single-cell RNA-seq by revisiting the first Alzheimer’s disease dataset. eLife 12, RP90214 (2023).

165. Hoffman, G. E. & Roussos, P. Dream: powerful differential expression analysis for repeated measures designs. Bioinformatics 37, 192–201 (2021).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

166. Wu, D. & Smyth, G. K. Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res. 40, e133 (2012).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

167. Viechtbauer, W. Conducting meta-analyses in R with the metafor package. J. Stat. Softw. 36, 1–48 (2010).

     Article  Google Scholar 

168. 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).

169. Bulik-Sullivan, B. K. et al. LD score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat. Genet. 47, 291–295 (2015).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

170. Chia, R. et al. Genome sequencing analysis identifies new loci associated with Lewy body dementia and provides insights into its genetic architecture. Nat. Genet. 53, 294–303 (2021).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

171. Ferrari, R. et al. Frontotemporal dementia and its subtypes: a genome-wide association study. Lancet Neurol. 13, 686–699 (2014).

     Article  PubMed  PubMed Central  Google Scholar 

172. Sealock, J. M. et al. Use of the PsycheMERGE network to investigate the association between depression polygenic scores and white blood cell count. JAMA Psychiatry 78, 1365–1374 (2021).

     Article  PubMed  PubMed Central  Google Scholar 

173. Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386 (2014).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

174. Bendall, S. C. et al. Single-cell trajectory detection uncovers progression and regulatory coordination in human B cell development. Cell 157, 714–725 (2014).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

175. Idrissi, B. Y., Arjovsky, M., Pezeshki, M. & Lopez-Paz, D. Simple data balancing achieves competitive worst-group-accuracy. Preprint at https://arxiv.org/abs/2110.14503 (2021).

176. Wang, J. Z., Du, Z., Payattakool, R., Yu, P. S. & Chen, C.-F. A new method to measure the semantic similarity of GO terms. Bioinformatics 23, 1274–1281 (2007).

     Article  CAS  PubMed  Google Scholar 

177. Zeng, B. et al. The single-cell and spatial transcriptional landscape of human gastrulation and early brain development. Cell Stem Cell 30, 851–866 (2023).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

178. Setty, M. et al. Characterization of cell fate probabilities in single-cell data with Palantir. Nat. Biotechnol. 37, 451–460 (2019).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

179. Donghoon, L. et al. Data for ‘Single-cell atlas of transcriptomic vulnerability across brain disorders’. Zenodo https://doi.org/10.5281/zenodo.14606776 (2025).

Download references

Author notes

1. These authors contributed equally: Donghoon Lee, Mikaela Koutrouli, Nicolas Y. Masse, Gabriel E. Hoffman

2. These authors jointly supervised this work: John F. Fullard, Panos Roussos

Authors and Affiliations

1. Center for Disease Neurogenomics, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Donghoon Lee, Mikaela Koutrouli, Nicolas Y. Masse, Gabriel E. Hoffman, Seon Kinrot, Xinyi Wang, Prashant N. M., Milos Pjanic, Tereza Clarence, Fotios Tsetsos, Deepika Mathur, David Burstein, Karen Therrien, Aram Hong, Clara Casey, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Jennifer Monteiro Fortes, Sarah R. Murphy, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, Hui Yang, Lyra Sheu, Marios Anyfantakis, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Sanan Venkatesh, Steven P. Kleopoulos, Zhenyi Wu, Kiran Girdhar, Georgios Voloudakis, Jaroslav Bendl, John F. Fullard & Panos Roussos

2. Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Donghoon Lee, Mikaela Koutrouli, Nicolas Y. Masse, Gabriel E. Hoffman, Seon Kinrot, Xinyi Wang, Prashant N. M., Milos Pjanic, Tereza Clarence, Fotios Tsetsos, Deepika Mathur, David Burstein, Karen Therrien, Aram Hong, Clara Casey, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Jennifer Monteiro Fortes, Sarah R. Murphy, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, Hui Yang, Lyra Sheu, Marios Anyfantakis, Maxim Signaevsky, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Sanan Venkatesh, Steven P. Kleopoulos, Zhenyi Wu, Kiran Girdhar, Georgios Voloudakis, Vahram Haroutunian, Jaroslav Bendl, John F. Fullard & Panos Roussos

3. Department of Psychiatry, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Donghoon Lee, Mikaela Koutrouli, Nicolas Y. Masse, Gabriel E. Hoffman, Seon Kinrot, Xinyi Wang, Prashant N. M., Milos Pjanic, Tereza Clarence, Fotios Tsetsos, Deepika Mathur, David Burstein, Karen Therrien, Aram Hong, Clara Casey, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Jennifer Monteiro Fortes, Sarah R. Murphy, Pavel Katsel, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, Hui Yang, Lyra Sheu, Marios Anyfantakis, Maxim Signaevsky, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Sanan Venkatesh, Steven P. Kleopoulos, Zhenyi Wu, Kiran Girdhar, Georgios Voloudakis, Vahram Haroutunian, Jaroslav Bendl, John F. Fullard & Panos Roussos

4. Department of Genetics and Genomic Sciences, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Donghoon Lee, Mikaela Koutrouli, Nicolas Y. Masse, Gabriel E. Hoffman, Seon Kinrot, Xinyi Wang, Prashant N. M., Milos Pjanic, Tereza Clarence, Fotios Tsetsos, Deepika Mathur, David Burstein, Karen Therrien, Aram Hong, Clara Casey, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Jennifer Monteiro Fortes, Sarah R. Murphy, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, Hui Yang, Lyra Sheu, Marios Anyfantakis, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Sanan Venkatesh, Steven P. Kleopoulos, Zhenyi Wu, Kiran Girdhar, Georgios Voloudakis, Jaroslav Bendl, John F. Fullard & Panos Roussos

5. Novo Nordisk Foundation Center for Protein Research, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark

   Mikaela Koutrouli, Lars J. Jensen & Lars Juhl Jensen

6. Center for Precision Medicine and Translational Therapeutics, James J. Peters VA Medical Center, New York, NY, USA

   Gabriel E. Hoffman, Fotios Tsetsos, Deepika Mathur, David Burstein, Sanan Venkatesh, Georgios Voloudakis & Panos Roussos

7. Mental Illness Research, Education and Clinical Center VISN2, James J. Peters VA Medical Center, New York, NY, USA

   Gabriel E. Hoffman, Deepika Mathur, David Burstein, Georgios Voloudakis, Vahram Haroutunian & Panos Roussos

8. Human Brain Collection Core, National Institute of Mental Health-Intramural Research Program, Bethesda, MD, USA

   Pavan K. Auluck & Stefano Marenco

9. Rush Alzheimer’s Disease Center, Rush University Medical Center, Chicago, IL, USA

   Lisa L. Barnes & David A. Bennett

10. Department of Neurological Sciences, Rush University Medical Center, Chicago, IL, USA

    Lisa L. Barnes & David A. Bennett

11. Department of Artificial Intelligence and Human Health, Icahn School of Medicine at Mount Sinai, New York, NY, USA

    Georgios Voloudakis

12. Department of Neuroscience, Icahn School of Medicine at Mount Sinai, New York, NY, USA

    Vahram Haroutunian

13. Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA

    Athan Z. Li, Daifeng Wang, Noah Cohen Kalafut & Sayali A. Alatkar

14. 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

15. Department of Biostatistics and Medical Informatics, University of Wisconsin-Madison, Madison, WI, USA

    Chenfeng He, Chirag Gupta, Daifeng Wang, Kalpana H. Arachchilage, Pramod B. Chandrashekar, Saniya Khullar & Ting Jin

16. Department of Psychiatry, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA

    Colleen A. McClung & Madeline R. Scott

17. Center for Systems and Therapeutics, Gladstone Institutes, San Francisco, CA, USA

    Gennadi Ryan, Monika Ahirwar, Steven Finkbeiner & Vivek G. Ramaswamy

18. Taube/Koret Center for Neurodegenerative Disease Research, Gladstone Institutes, San Francisco, CA, USA

    Gennadi Ryan, Monika Ahirwar, Steven Finkbeiner & Vivek G. Ramaswamy

19. Department of Population Health Sciences, University of Wisconsin-Madison, Madison, WI, USA

    Jerome J. Choi

20. Vanderbilt Genetics Institute, Vanderbilt University Medical Center, Nashville, TN, USA

    Logan C. Dumitrescu & Timothy J. Hohman

21. Vanderbilt Memory & Alzheimer’s Center, Vanderbilt University Medical Center, Nashville, TN, USA

    Logan C. Dumitrescu & Timothy J. Hohman

22. Department of Neurology, University of California San Francisco, San Francisco, CA, USA

    Steven Finkbeiner

23. Department of Physiology, University of California San Francisco, San Francisco, CA, USA

    Steven Finkbeiner

24. Neuroscience and Biomedical Sciences Graduate Programs, University of California San Francisco, San Francisco, CA, USA

    Steven Finkbeiner

Consortia

Contributions

Conceptualization: D.L., V.H., J.F.F. and P.R. Methodology: D.L., G.E.H., J.F.F. and P.R. Software: D.L., N.Y.M. and G.E.H. Validation: D.L., X.W., J.M.F., S.R.M., P.K. and J.F.F. Formal analysis: D.L., M.K., N.Y.M., S.K., M.P., T.C., F.T. and G.E.H. Investigation: X.W., A.H., C.C., Z.S., M.A., S.A., J.M.F. and J.F.F. Resources: P.K.A., L.L.B., S.M., D.A.B. and V.H. Data curation: D.L., P.N.M., D.M., D.B., K.T. and J.B. Writing: D.L., N.Y.M., J.F.F. and P.R., with support from all of the authors. Visualization: D.L., M.K. and N.Y.M. Supervision: D.L., L.J.J., K.G., G.V., V.H., J.B., J.F.F. and P.R. Project administration: D.L. and P.R. Funding acquisition: D.L., V.H. and P.R. All of the authors read and approved the final draft of the paper.

Corresponding authors

Correspondence to Donghoon Lee or Panos Roussos.

Competing interests

The authors declare no competing interests.

Peer review information

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

Open Access This article is licensed under a Creative Commons Attribution 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