Polygenic enrichment distinguishes disease associations of individual cells in single-cell RNA-seq data

Single-cell RNA sequencing (scRNA-seq) provides unique insights into the pathology and cellular origin of disease. We introduce single-cell disease relevance score (scDRS), an approach that links scRNA-seq with polygenic disease risk at single-cell resolution, independent of annotated cell types. scDRS identifies cells exhibiting excess expression across disease-associated genes implicated by genome-wide association studies (GWASs). We applied scDRS to 74 diseases/traits and 1.3 million single-cell gene-expression profiles across 31 tissues/organs. Cell-type-level results broadly recapitulated known cell-type–disease associations. Individual-cell-level results identified subpopulations of disease-associated cells not captured by existing cell-type labels, including T cell subpopulations associated with inflammatory bowel disease, partially characterized by their effector-like states; neuron subpopulations associated with schizophrenia, partially characterized by their spatial locations; and hepatocyte subpopulations associated with triglyceride levels, partially characterized by their higher ploidy levels. Genes whose expression was correlated with the scDRS score across cells (reflecting coexpression with GWAS disease-associated genes) were strongly enriched for gold-standard drug target and Mendelian disease genes. scDRS associates individual cells in scRNA-seq with disease by scoring single-cell transcriptomes using GWAS gene signatures. Applied to 74 GWAS and 1.3 million single-cell profiles, scDRS identifies specific cellular subpopulations associated with these diseases.

T he mechanisms through which risk variants identified by GWASs impact critical tissues and cell types are largely unknown 1 ; identifying these tissues and cell types is central to our understanding of disease etiologies and will inform efforts to develop therapeutic treatments 2 . scRNA-seq has emerged as the tool of choice for measuring gene expression abundance at single-cell resolution 3 , providing an increasingly clear picture of which genes are active in which cell types and identifying new cell populations within classically defined cell types. Integrating scRNA-seq with GWAS data offers the potential to identify critical tissues, cell types and cell populations through which GWAS risk variants impact disease 4-6 , thus providing finer resolution than studies using bulk transcriptomic data [7][8][9][10] .
Previous studies integrating scRNA-seq with GWAS have largely focused on predefined cell-type annotations (for example, classical cell types defined using known marker genes), aggregating cells from the same cell type followed by evaluating overlap of the cell-type-level information with GWAS 4-6 . However, this approach overlooks the considerable heterogeneity of individual cells within cell types that has been reported in studies of scRNA-seq data alone [11][12][13][14][15][16] ; the underlying methods (for example, Seurat cell-scoring function 13 , Vision 14 and VAM 16 ) have sought to explain transcriptional heterogeneity in scRNA-seq data by scoring cells based on predefined gene sets such as cellular pathways but do not consider polygenic disease risk from GWASs and generally do not provide individual-cell-level association P values. Integrating information from scRNA-seq data at fine-grained resolution (for example, individual cells both within and across cell types) with polygenic signals from disease GWAS has considerable potential to produce new biological insights.
Here, we introduce the single-cell disease relevance score (scDRS), a method to evaluate polygenic disease enrichment of individual cells in scRNA-seq data. scDRS assesses whether a given cell has excess expression levels across a set of putative disease genes derived from GWASs using an appropriately matched empirical null distribution to estimate well-calibrated P values. We performed extensive simulations to assess the calibration and power of scDRS. We applied scDRS to 74 diseases/traits (average GWAS N = 346K where "K" means thousands) and 16 scRNA-seq data sets (including the Tabula Muris Senis (TMS) mouse cell atlas 17 ), assessing cell-typedisease associations and within-cell-type association heterogeneity, including heterogeneity of T cells in association with autoimmune diseases, neurons in association with brain-related diseases/traits and hepatocytes in association with metabolic traits. We analyzed a broader set of scRNA-seq data sets to provide validation across species (human versus mouse) and across sequencing platforms and include scRNA-seq data sets with experimentally determined cell types and cell states.

Results
Overview of methods. scDRS integrates gene expression profiles from scRNA-seq with polygenic disease information from GWASs to associate individual cells to disease, by assessing the excess expression of GWAS putative disease genes in a given cell relative to other genes with similar expression across all cells. scDRS consists of three steps (Fig. 1). First, scDRS constructs a set of putative disease genes from GWAS summary statistics using MAGMA 18 , an existing gene scoring method (top 1,000 MAGMA genes). Second, scDRS quantifies the aggregate expression of the putative disease genes in each cell to generate cell-specific raw disease scores; to maximize power, each putative disease gene is weighted by its GWAS MAGMA z-score and inversely weighted by its gene-specific technical noise level in the single-cell data, estimated via modeling the mean-variance relationship across genes 16,19 . To determine statistical significance, scDRS also generates 1,000 sets of cell-specific raw control scores at Monte Carlo (MC) samples of matched control gene sets (matching gene set size, mean expression and expression variance of the putative disease genes). Third, scDRS normalizes the raw disease score and raw control scores for each cell (producing the normalized disease score and normalized control scores), and then computes cell-level P values based on the empirical distribution of the pooled normalized control scores across all control gene sets and all cells. Further details are provided in Methods, Supplementary Note and Supplementary Figs. 1-3. scDRS outputs individual-cell-level disease-association P values, normalized disease scores, and 1,000 sets of normalized control scores (referred to as "disease scores" and "control scores" in the rest of the paper) that can be used for a wide range of downstream applications (Methods). Here, we focus on three downstream analyses. First, we perform cell-type-level analyses to associate predefined cell types with disease and assess heterogeneity in association with disease across cells within a predefined cell type. Second, we perform individual-cell-level analyses to associate individual cells to disease and correlate individual-cell-level variables to the scDRS (1) scDRS constructs a set of putative disease genes from GWAS summary statistics by selecting the top 1,000 MAGMA genes; these putative disease genes are expected to have higher expression levels in the relevant cell population. (2) scDRS computes a raw disease score for each cell, quantifying the aggregate expression of the putative disease genes in that cell; to maximize power, each putative disease gene is weighted by its GWAS MAGMA z-score and inversely weighted by its gene-specific technical noise level in scRNA-seq. scDRS also computes a set of 1,000 MC raw control scores for each cell, in each case using a random set of control genes matching the gene set size, mean expression and expression variance of the putative disease genes. (3) scDRS normalizes the raw disease score and raw control scores across gene sets and across cells and then computes a P value for each cell based on the empirical distribution of the pooled normalized control scores across all control gene sets and all cells. The choice of 1,000 for the number of putative disease genes and the choice of 1,000 for the number of control scores are independent.
disease score. Third, we perform gene-level analyses to prioritize disease-relevant genes whose expression is correlated with the scDRS disease score, reflecting coexpression with genes implicated by disease GWASs. We analyzed publicly available GWAS summary statistics of 74 diseases/traits (average N = 346K; Supplementary Table 1) in conjunction with 16 scRNA-seq or single-nucleus RNA sequencing (snRNA-seq) data sets spanning 1.3 million cells from 31 tissues/organs from mouse (Mus musculus) and human (Homo sapiens) (Supplementary Tables 2-7 and Supplementary Fig. 4). The single-cell data sets include two data sets from the TMS mouse cell atlases 17 collected using different technologies, the Tabula Sapiens (TS) human cell atlas 20 , and other data sets focusing on specific tissues containing well-annotated cell types/states. We focused on the TMS fluorescence-activated cell sorting (FACS) data in our primary analyses due to its comprehensive coverage of 23 tissues and 120 cell types and more accurate quantification of gene expression levels (via SMART-Seq2 21 ); we used the other 15 data sets to validate our results. We note the extensive use of mouse gene expression data to study human diseases and complex traits 4-7,10,22 (Supplementary Note).
Simulations assessing calibration and power. We performed null simulations and causal simulations to assess the calibration and power of scDRS, comparing scDRS to three state-of-art methods for scoring individual cells with respect to a specific gene set: Seurat (cell-scoring function) 13 , Vision 14 and VAM 16 (Methods).
First, we evaluated each method in null simulations where no cells have systematically higher expression across the putative disease genes analyzed. We subsampled 10,000 cells from the TMS FACS data and randomly selected 1,000 putative disease genes. We simulated GWAS gene weights for scDRS matching the MAGMA z-score distributions in real traits and used binary disease gene sets for the other three methods. scDRS and Seurat produced well-calibrated P values, Vision showed slightly inflated type I errors and VAM produced severely inflated type I errors (Fig. 2a  and Supplementary Table 8).
Next, we evaluated scDRS, Seurat and Vision in causal simulations where a subset of causal cells has systematically higher expression across putative disease genes (we did not include VAM, which was not well calibrated in null simulations). We used the same 10,000 cells subsampled from the TMS FACS data, randomly selected 1,000 causal disease genes, randomly selected 500 of the 10,000 cells as causal cells and artificially perturbed their expression levels to be higher (1.05-1.50× for different simulations) across the 1,000 causal disease genes, and we randomly selected 1,000 putative disease genes (provided as input to each method) with 25% overlap with the 1,000 causal disease genes. We used the binary gene set for all three methods because there were no GWAS weights involved in data generation. We determined that scDRS attained higher power than Seurat and Vision to detect individual-cell-disease associations at false discovery rate (FDR) < 0.1 ( Fig. 2b and Supplementary  Table 9); the improved power of scDRS may be due to its incorporation of gene-specific weights that downweight genes with higher levels of technical noise.
Results of additional null and causal simulations are reported in the Supplementary Note, Extended Data Figs. 1, 2, Supplementary  Fig. 5 and Supplementary Table 10.
Results across 120 TMS cell types for 74 diseases and traits. We analyzed GWAS data from 74 diseases/traits (average N = 346K; Supplementary Tables 1 and 11) in conjunction with the TMS FACS data with 120 cell types (cells from different tissues were combined for a given cell type; Supplementary Table 5). We first report scDRS cell-type-level results, aggregated for each cell type from the scDRS individual-cell-level results; the individual-cell-level results are discussed in subsequent sections. Results for a representative subset of 19 cell types and 22 diseases/traits are reported in Fig. 3. Within this subset, scDRS identified 80 associated cell-typedisease pairs (FDR < 0.05; squares in Fig. 3) and detected significant within-cell-type disease-association heterogeneity for 43 of these 80 associated cell-type-disease pairs (FDR < 0.05; cross symbols in Fig. 3) (273 of 597 across all pairs of the 120 cell types and  74 diseases/traits; Extended Data Figure 3 and Supplementary  Table 12).
For cell-type-disease associations, as expected, scDRS broadly associated blood/immune cell types with blood/immune-related diseases/traits, brain cell types with brain-related diseases/traits and other cell types with other diseases/traits (block-diagonal pattern in Fig. 3). Most scDRS discoveries recapitulated well-established cell-type-disease associations, including blood/immune cell types with blood cell traits, immune cell types with immune diseases, neuronal cell types with brain-related traits/diseases 6,22,23 and hepatocytes with metabolic traits 24 . In addition, chondrocytes, bladder cells, ventricular myocytes and pancreatic beta cells were associated with their corresponding expected diseases/traits [25][26][27][28] . Exceptions to the block-diagonal pattern and further details are discussed in the Supplementary Note.
We also discuss several less well-established results. First, granulocyte monocyte progenitors were strongly associated with multiple sclerosis (MS), highlighting the role of myeloid cells in MS 29 . Second, oligodendrocytes and oligodendrocyte precursor cells (OPCs) were associated with multiple brain-related diseases/traits. These associations are less clear in existing genetic studies 4,6,22,30 but are biologically plausible, consistent with the increasingly discussed role of oligodendrocyte lineage cells in brain diseases/traits; the differentiation and myelination activity of oligodendrocyte lineage cells are important to maintain the functionality of neuronal cells 31 . We detected significant heterogeneity across OPCs in their association with many brain-related diseases/traits, consistent with recent evidence of functionally diverse states of OPCs 32 , traditionally considered to be a homogeneous population ( Supplementary  Fig. 6). Third, we detected significant heterogeneity across cells for the association between proerythroblasts and red cell distribution width, which may correspond to erythrocytes at different differentiation stages 33 (Supplementary Fig. 7). We also detected other instances of significant within-cell-type association heterogeneity, including T cells with immune diseases, neurons with brain-related diseases/traits and hepatocytes with metabolic traits (discussed in subsequent sections). Additional

Heterogeneous T cells subpopulations in autoimmune disease.
We investigated the heterogeneity across TMS FACS T cells in association with autoimmune diseases (Fig. 3). We jointly analyzed all TMS FACS T cells (3,769 cells spanning 15 tissues). Because the original study clustered cells from different tissues separately 17 , we reclustered these T cells, resulting in 11 clusters ( Fig. 4a and Methods); we verified that batch effects were not observed for tissue, age or sex ( Supplementary Fig. 12). We considered ten autoimmune diseases: inflammatory bowel disease (IBD), Crohn's disease (CD), ulcerative colitis (UC), rheumatoid arthritis (RA), MS, autoimmune traits (AITs; a general term for autoimmune diseases), hypothyroidism (HT), eczema, asthma (ASM) and respiratory and ear-nose-throat diseases (RR-ENT) (Supplementary Table 1); we considered height as a negative control trait.
We focused on individual cells associated with IBD, a representative autoimmune disease (Fig. 4b). The 387 IBD-associated cells (FDR < 0.1) formed subpopulations of 4 of the 11 T cell clusters. The subpopulation of 123 IBD-associated cells in cluster 3 (labeled T reg -like) had high expression of regulatory T (T reg ) cell marker genes (FOXP3 + , CTLA4 + and LAG3 + ; Supplementary Fig. 17a), and their specifically expressed genes significantly overlapped with T reg signatures (P = 6.0 × 10 −8 for MSigDB signatures and P = 4.0 × 10 −68 for an effector-like T reg program 34 , two-sided Fisher's exact test; Supplementary Fig. 17c,d), suggesting these cells had T reg immunosuppressive functions. Interestingly, these 123 IBD-associated cells were nonrandomly distributed in cluster 3 on the uniform manifold approximation and projection (UMAP) plot (P < 0.001, MC test; Methods). Genes specifically expressed in these IBD-associated cells were preferentially enriched (compared to the 506 non-IBD-associated cells in the same cluster) in pathways defined by nuclear factor κB signaling, T helper cell differentiation and tumor necrosis factor-mediated signaling ( Supplementary Fig. 17e); these pathways are closely related to inflammation, a distinguishing feature of IBD 35 . In addition, the subpopulations of IBD-associated cells in clusters 4, 5 and 9 were labeled as "Th2/T reg -like", "Th17-like" and "CD8 + effector-like", respectively, consistent with previous studies associating subpopulations of effector T cells to IBD, particularly T reg cells and Th17 cells 35 . Results for other autoimmune diseases, details for annotating disease-associated T cell subpopulations and replication analyses on two human scRNA-seq data sets 36,37 are reported in Fig. 4c We further compared the individual T cell associations of IBD with HT, another representative autoimmune disease (Fig. 4c,d). The top four HT-associated subpopulations included not only three IBD-associated subpopulations (cells in clusters 3, 4 and 9; Fig. 4c) but also a unique subpopulation of cells in cluster 10 (labeled as "proliferative"). Despite the stronger associations with HT overall (possibly due to higher GWAS power), IBD was more strongly associated with cells in cluster 4 (labeled as "Th2/T reg -like";  Table 19. Motivated by the effector-like T cell subpopulations associated with IBD, we investigated whether the heterogeneity of T cells in association with autoimmune diseases was correlated with T cell effectorness gradient, a continuous classification from naive to effector T cells. We separately computed the effectorness gradients for CD4 + T cells (1,686 cells) and CD8 + T cells (2,197 cells) using pseudotime analysis 36,38 and assessed whether the CD4 (respectively, CD8) effectorness gradient was correlated with scDRS disease scores for the ten autoimmune diseases across CD4 + T cells (respectively, CD8 + T cells). Results are reported in Fig. 4e and Supplementary Table 20. The CD4 effectorness gradient had strong associations with IBD, CD, UC, AIT and HT (P < 0.005, MC test), weak associations with eczema and ASM (P < 0.05) and nonsignificant associations with RA, MS and RR-ENT, implying these autoimmune diseases are associated with more effector-like CD4 + T cells. The CD8 effectorness gradient had weaker associations (P < 0.05 for IBD, CD and AIT but nonsignificant for others), suggesting that CD4 + effector T cells may be more important than CD8 + effector T cells for these diseases. The association of T cell effectorness gradients with autoimmune diseases has not been formally evaluated previously but is consistent with previous studies linking T cell effector functions to autoimmune disease 39 or characterizing similarities among effector T cell subtypes 36,40 . Additional results on T cell effectorness gradient and comparison to cluster-level LDSC-SEG 10 Tables 17 and 20.
Finally, we prioritized disease-relevant genes by computing the correlation (across all TMS FACS cells) between the expression of a given gene and the scDRS score for a given disease; this approach identifies genes coexpressed with genes implicated by disease GWAS. We compared the top 1,000 genes prioritized using this approach with putative drug targets 41 (for eight autoimmune diseases and excluding RR-ENT and HT) or genes known to cause a Mendelian form of the disease 42   reported in Fig. 4f and Supplementary Table 22. We determined that scDRS attained a more accurate prioritization of disease-relevant genes compared to the top 1,000 MAGMA genes (2.07 for median ratio of (excess overlap − 1); Methods), likely by capturing disease-relevant genes with weak GWAS signal 43 . For example, scDRS prioritized ITGB7 for IBD (rank 11; drug target for IBD using vedolizumab 44 ) and JAK1 for RA (rank 358; drug target for RA using baricitinib 45 ), both missed by MAGMA (ranks 10,565 and 5,228, P = 0.54, 0.26, respectively). Additional results are reported in the Supplementary Note, Extended Data Figure 7 and Supplementary  Table 21.
Heterogeneous neuron subpopulations in brain traits. We investigated the heterogeneity across neurons in association with brain-related diseases/traits (Fig. 3). We considered seven brainrelated diseases/traits: schizophrenia (SCZ), major depressive disorder (MDD), bipolar disorder (BP), neuroticism (NRT), college education (ECOL), body mass index (BMI) and smoking (Supplementary Table 1); we considered height as a negative control trait. Because the TMS FACS data had limited coverage of neuronal subtypes, we focused on a separate mouse brain scRNA-seq data (Zeisel & Muñoz-Manchado et al. 46 ; 3,005 cells), which had better coverage of neuronal subtypes and have been analyzed at the cell-type Systolic blood pressure (SBP)   Table 1). Abbreviated cell-type names include red blood cells, granulocyte monocyte progenitors, medium spiny neuron (MSN) and OPCs. Neuron refers to neuronal cells with undetermined subtypes (whereas MSN and interneuron (nonoverlapping with neuron) refer to neuronal cells with those inferred subtypes). Complete results for 120 cell types and 74 diseases/traits are reported in Extended Data  . e, Association between scDRS disease score and CD4 effectorness gradient across CD4 + T cells for five representative autoimmune diseases and height, a negative control trait. x axis denotes CD4 effectorness gradient quintile bins, and y axis denotes the average scDRS disease score in each bin for each disease. *P < 0.05 and, **P < 0.005 (one-sided MC test). Numerical results for all ten autoimmune diseases are reported in Supplementary Table 20. f, Excess overlap with gold-standard gene sets. x axis denotes excess overlap of genes prioritized by MAGMA, and y axis denotes excess overlap of genes prioritized by scDRS, for each of the ten autoimmune diseases. Median ratio of (excess overlap − 1) for scDRS versus MAGMA is 2.07. Numerical results are reported in Supplementary Table 22.  Fig. 5a (results for all seven brain-related traits are shown in Supplementary Fig. 23b). We observed a continuous gradient of CA1 pyramidal neuron associations with the seven brain-related traits. Motivated by known location-specific functions of hippocampal neurons 15 , we investigated whether the heterogeneity observed in Fig. 5a was correlated with spatial location. We considered the natural CA1 spatial axes 48 (dorsal-ventral long axis, proximal-distal transverse axis and superficial-deep radial axis) and inferred spatial coordinates for each cell in terms of six continuous individual-cell-level scores for these spatial regions ( Supplementary  Figs. 23c and 24, Supplementary Table 24 and Methods). The inferred spatial scores for the dorsal-ventral and proximal-distal axes varied along the top two UMAP axes, providing visual evidence of stronger neuron-SCZ associations in dorsal and proximal regions ( Fig. 5a and Supplementary Fig. 23).

Supplementary
We used the results of scDRS for individual cells to assess whether the inferred spatial scores for each of the six spatial regions (dorsal/ ventral/proximal/distal/superficial/deep) were correlated to the scDRS disease scores for each of the seven brain-related traits across CA1 pyramidal neurons (Methods). Results are reported in Fig. 5b scDRS SCZ score (CA1 pyramidal neurons) ( . We include a visualization of putative dorsal-ventral and proximal-distal axes (see the "Heterogeneous neuron subpopulations in brain traits" subsection). Results for all seven brain-related diseases/traits and height are reported in Supplementary Fig. 23b. b, Association between scDRS disease score and proximal score across CA1 pyramidal neurons for five representative brain-related disease/traits and height, a negative control trait. x axis denotes proximal score quintile bins, and y axis denotes average scDRS disease score in each bin for each disease. *P < 0.05, **P < 0.005 (one-sided MC test). Results for all six spatial scores and all seven brain traits (and height) are reported in Extended Data Figure 8 and Supplementary Table 25. BP, bipolar disorder; ECOL, college education; MDD, major depressive disorder. c, Subpopulations of hepatocytes associated with triglyceride levels (TG) in the TMS FACS data. Significantly associated cells (FDR < 0.1) are denoted in red, with shades of red denoting scDRS disease scores; nonsignificant cells are denoted in gray. Cluster boundaries indicate the corresponding hepatocyte clusters. In the legend, numbers in parentheses denote the number of TG-associated cells versus the total number of cells. Cluster labels are based on the putative identity of cells in the cluster. Results for the other eight metabolic traits and height are reported in Supplementary Fig. 25. d, Association between scDRS disease score and polyploidy score for four representative metabolic traits and height, a negative control trait. x axis denotes polyploidy score quintile bins, and y axis denotes average scDRS disease score in each bin for each disease. *P < 0.05, **P < 0.005 (one-sided MC test). Results for all three scores (polyploidy score, pericentral score and periportal score) and all nine metabolic traits (and height) are reported in Extended Data Figure 9 and Supplementary Table 26. HDL, high-density lipoprotein; TC, total cholesterol; TST, testosterone.
(for the proximal region, which had the strongest associations), Extended Data Figure 8 and Supplementary Table 25. The proximal score was strongly associated with all seven brain-related traits (all P < 0.002, MC test), suggesting proximal CA1 pyramidal neurons may be more relevant to these brain-related traits (instead of distal CA1 pyramidal neurons). The association between the proximal region and brain-related traits is consistent with the fact that the proximal region of the hippocampus receives synaptic inputs in the perforant pathway, which is the main input source of the hippocampus 49 . Validations of the spatial scores using independent data and results on other spatial scores and six additional mouse and human data sets [50][51][52][53][54][55] are reported in the Supplementary Note and Extended Data Figure 8.

Heterogeneous hepatocyte subpopulations in metabolic traits.
We investigated the heterogeneity across TMS FACS hepatocytes (in the liver) in their association with metabolic traits (Fig. 3). We considered nine metabolic traits: total cholesterol (TG), high-density lipoprotein (HDL), low-density lipoprotein (LDL), total cholesterol (TC), testosterone (TST), alanine aminotransferase (ALT), alkaline phosphatase (ALP), sex hormone-binding globulin (SHBG) and total bilirubin (TBIL) (Supplementary Table 1); we considered height as a negative control trait. We focused on individual cells associated with TG, a representative metabolic trait ( Fig. 5c; results for other traits in Supplementary  Fig. 25). The 530 TG-associated cells (FDR < 0.1) formed subpopulations of five of the six hepatocyte clusters; we characterized these subpopulations based on ploidy level (number of sets of chromosomes in a cell) and zonation (pericentral/mid-lobule/periportal spatial location in the liver lobule), which have been extensively investigated in previous studies of hepatocyte heterogeneity 56,57 . We inferred the ploidy level and zonation for each individual cell in terms of continuous individual-cell-level polyploidy, pericentral and periportal scores (Supplementary Fig. 26 and Methods). The inferred ploidy level and zonation varied across clusters, providing visual evidence of stronger cell-TG associations in high-ploidy clusters (clusters 1 and 2), particularly the periportal high-ploidy cluster (cluster 2; Fig. 5c).
We used the results of scDRS for individual cells to assess whether the inferred polyploidy, pericenteral and periportal scores were correlated to the scDRS disease score for each of the nine metabolic traits across hepatocytes; we jointly regressed the scDRS disease score for each trait on the polyploidy score, pericentral score and periportal score (because the polyploidy score was positively correlated with the other two scores; Methods). Results are reported in Fig. 5d (for the polyploidy score which had the strongest associations), Extended Data Figure 9 and Supplementary Table 26. The polyploidy score was strongly associated with all nine metabolic traits (all P < 0.007, MC test), suggesting that high-ploidy hepatocytes may be more relevant to these metabolic traits. The association between ploidy level and metabolic traits is consistent with previous findings that ploidy levels are associated with changes in the expression level of genes for metabolic processes such as de novo lipid biosynthesis and glycolysis 57,58 and supports the hypothesis that liver functions are enhanced in polyploid hepatocytes 57 . In addition, the periportal score was associated with the nine metabolic traits (all P < 0.005 except P = 0.04,0.02,0.24 for HDL, SHBG and TBIL, MC test). Although the pericentral score was not significantly associated with these traits in TMS FACS, we detected significant associations across multiple other data sets (Supplementary Note). These results suggest that these metabolic traits are impacted by complex processes involving both pericentral and periportal hepatocytes. Validations of the polyploidy and zonation scores using independent data and results on five additional mouse and human data sets 17

Discussion
We have introduced scDRS, a method that leverages polygenic GWAS signals to associate individual cells in scRNA-seq data with diseases and complex traits. We showed via extensive simulations that scDRS is well calibrated and powerful in realistic scenarios. We applied scDRS to 74 diseases/traits in conjunction with 16 scRNA-seq data sets and detected extensive heterogeneity in disease associations of individual cells within classical cell types. These findings may prove useful for targeting the relevant cell populations for in vitro experiments to elucidate the molecular mechanisms through which GWAS risk variants impact disease.
We have demonstrated the value in associating individual cells to disease, assessing the heterogeneity across individual cells within predefined cell types in their association to disease, identifying cell-level variables partially characterizing the individual cells that are associated to disease, and broadly associating predefined cell types to disease. Analyses of larger scRNA-seq/snRNA-seq and GWAS data sets using approaches such as scDRS will continue to further these goals.
We note several limitations and future directions of our work. First, identifying a statistical correlation between individual cells (or cell types) and disease does not imply causality but may instead reflect indirect tagging of causal cells/cell types, analogous to previous works 4,5,10,18 . However, even in such cases, the implicated cells/ cell types are likely to be closely biologically related to the causal cells/cell types, based on their similar expression patterns. Second, the fact that scDRS assesses the statistical significance of an individual cell's association to disease by implicitly comparing it to other cells via matched control genes may reduce power if most cells in the data are truly causal. For example, association with IBD in a data set containing only T reg cells (one of the causal cell types for IBD) will likely yield largely nonsignificant results. This limitation did not impact our main analyses, because the TMS data include a broad set of cell types; in more specialized data sets (which may be preferred in some settings due to the more comprehensive profiling of the focal cell population), this limitation can potentially be addressed by selecting matched control genes based on a broad cell atlas (for example, the TMS or TS data). Third, although we have primarily focused on the associations involving a single disease/trait, further investigation of differences between diseases/traits within the same category is an important future direction. Additional limitations are discussed in the Supplementary Note. Despite all these limitations, scDRS is a powerful method for distinguishing disease associations of individual cells in scRNA-seq data.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41588-022-01167-z.

Methods
Ethics statement. This study analyzed publicly available data sets and hence did not require ethical approval. scDRS method. We consider an scRNA-seq data set with n cell cells (not cell types) and n gene genes. We denote the cell-gene matrix as X ∈ R ncell×ngene , where X cg represents the expression level of cell c and gene g. We assume that X is size-factor-normalized (for example, 10,000 counts per cell) and log-transformed (log(x + 1)) from the original raw count matrix 19 . We regress the covariates out from the normalized data 19 (with a constant term in the regressors to center the data), before adding the original log mean expression of each gene back to the residual data. Such a procedure preserves the mean-variance relationship in the covariate-corrected data, necessary for estimating the gene-specific technical noise levels (Supplementary Note). Gene-level statistics for the scRNA-seq data are reported in Supplementary  Fig. 4 and Supplementary Tables 3 and 4; estimated technical noise levels are moderately correlated across genes between the 16 data sets (average correlation = 0.34) and are highly correlated between data sets with similar cell-type compositions (for example, 0.74 between TMS FACS and TS FACS).
The scDRS algorithm is described below. Given a disease GWAS and an scRNA-seq data set, scDRS computes a P value for each individual cell for association with the disease. scDRS also outputs cell-level normalized disease scores and B sets of normalized control scores (default B = 1,000) that can be used for data visualization and MC-based statistical inference (see below). scDRS consists of three steps. First, scDRS constructs a set of putative disease genes from the GWAS summary statistics. Second, scDRS computes a raw disease score and B MC samples of raw control scores for each cell. Third, after gene set-wise and cell-wise normalization, scDRS computes an association P value for each cell by comparing its normalized disease score to the empirical distribution of the pooled normalized control scores across all control gene sets and all cells.
We discuss guidelines for using scDRS. First, scDRS relies on assumptions that require the disease gene set to have a moderate size (for example, >50 genes and <20% of all genes). Second, to ensure a reasonable number of scDRS discoveries, we recommend using GWAS data with a heritability z-score greater than 5 or a sample size greater than 100K. We also recommend using scRNA-seq data with a diverse set of cells potentially relevant to disease, although a smaller number of cells should not affect the scDRS power. However, scDRS will not produce false positives for less ideal GWAS or single-cell data sets. Third, scDRS is computationally efficient and scales linearly with the number of cells and number of control gene sets for both computation and memory use; it takes around 3 h and 60 GB for a single-cell data set with a million cells). Fourth, scDRS can be used in conjunction with any gene sets (instead of the MAGMA GWAS gene sets) to evaluate the enrichment of aggregate expression for cells in a single-cell data set. Fifth, we provide an option to adjust for cell-type proportions (or any cell group annotations), recommended only for extremely unbalanced data sets. Further details and evaluations of alternative versions of scDRS are provided in the Supplementary Note.

Compute disease-association P values
a. First gene set alignment by mean and variance. Let σ 2 g be the expression variance of gene g. For each cell c, b. Cell-wise standardization for each cell c by the mean μ ctrl c and SD σ ctrl c of control scores s ctrl c1 , …, s ctrl cB of that cell, c. Second gene set alignment by mean. For each cell c, Compute cell-level P values based on the empirical distribution of the pooled normalized control scores for each cell c, Output: cell-level P values p c , normalized disease scores s c , and normalized control scores s ctrl c1 , …, s ctrl cB .

Downstream applications and MC test.
We use a unified MC test for the scDRS downstream analyses based on the (normalized) disease and control scores. Specifically, let t be a test statistic computed from the disease score of the given set of cells (different downstream analyses differ by the test statistics), and let t ctrl 1 , …, t ctrl B be the same test statistics computed from the B sets of control scores of the same set of cells. The MC P value can be written as: The MC test avoids the assumption that the cells are independent-a strong assumption in scRNA-seq analyses, for example, when analyzing cells in the same cluster that are dependent due to the clustering process. However, the MC P value p MC cannot be smaller than 1/(1 + B) by equation (6), limiting its ability in distinguishing highly significant associations. We can also compute an MC z-score as z MC ; this MC z-score is not restricted by the MC limit of 1/(1 + B) but relies on the assumption that the control test statistics approximately follow a normal distribution. We recommend using MC P values to determine statistical significance and using MC z-scores to further prioritize associations whose MC P values have reached the MC limit. Below, we describe the test statistics used by the three analyses listed above. Besides the three analyses below, the MC test can in principle be extended to any analysis that computes a test statistic from the disease score of a set of cells.
Associating cell type with disease. We use the top 5% quantile of the disease score of cells from the given cell type as the test statistic. This test statistic is robust to annotation outliers, for example, a few misannotated but highly significant cells. One can also use other test statistics such as the top 1% quantile or the maximum.
Assessing within-cell-type heterogeneity in association with disease. We use Geary's C 14,63 as the test statistic. Geary's C measures the spatial autocorrelation of the disease score across a set of cells (for example, cells from the same cell type or cell cluster) with respect to a cell-cell similarity matrix. Given a set of n cells, the corresponding disease scores s 1 ,…s n , and the cell-cell similarity matrix W ∈ R n×n , Geary's C is calculated as where s = 1 n ∑ n i=1 si. A value significantly lower than 1 suggests a high level of disease-association heterogeneity across the given set of cells. Details are provided in the Supplementary Note.

Correlating a cell-level variable to disease across a given set of cells.
For associating a single-cell-level variable with disease, we use the Pearson's correlation between the cell-level variable and the disease score across the given set of cells as the test statistic. For jointly associating multiple cell-level variables with disease, we use the regression t-statistic as the test statistic, obtained from jointly regressing the disease score against the cell-level variables.
Simulations. We considered three comparison methods: Seurat 13 ("score_genes" as implemented in Scanpy 64 v1.6.0), Vision 14 (implemented in scDRS v1.0.1) and VAM 16 (v0.5.1). To our knowledge, VAM is the only published cell-scoring method that provides cell-level association P values. We chose to include Seurat due to its wide use. We standardized its output cell-level scores (mean 0 and SD 1) before computing the cell-level P values based on the standard normal distribution.
We chose to include Vision because its outputs are nominal cell-level z-scores that can be easily converted to P values; we again added the standardization step, because otherwise, the Vision results were highly unstable. We did not include other methods like PAGODA 11 or AUCell 12 , because it is not straightforward to convert their outputs to cell-level association P values and also because the z-scoring methods (for example, Vision) outperformed other methods in a recent evaluation 16 .
We performed simulations on a data set with 10,000 cells subsampled from the TMS FACS data. In null simulations, we randomly selected putative disease genes from a set of noninformative genes. We considered four numbers of putative disease genes (100, 500, 1,000 or 2,000) and four types of genes to sample from: (1) the set of all genes, (2) the set of top 25% genes with high mean expression, (3) the set of top 25% genes with high expression variance and (4) the set of top 25% overdispersed genes, where the level of overdispersion is calculated as the difference between the actual variance and the estimated technical variance in the log scale data. For the default version of scDRS, we simulated GWAS gene weights by first randomly selecting a disease (out of the 74 diseases/traits) and then randomly permuting the top MAGMA z-scores from the selected disease. We did not simulate gene-specific technical noise-based single-cell weights, because these weights were inherent to the single-cell data. For the MC test for cell-type-disease association, we used the top 5% quantile as the test statistic (see above).
In causal simulations, we randomly selected 1,000 causal disease genes, randomly selected 500 of the 10,000 cells as causal cells and artificially perturbed their expression levels to be higher (at various effect sizes) across the 1,000 causal disease genes and randomly selected 1,000 putative disease genes (provided as input to scDRS and other methods) with various levels of overlap with the 1,000 causal disease genes. Here, the effect size corresponds to the fold change of expression of the causal genes in the causal cells (multiplicative in the original count space and additive in the log space). We performed three sets of causal simulations: (1) varying effect size from 5% to 50% while fixing 25% overlap, (2) varying level of overlap from 5% to 50% while fixing 25% effect size and (3) assigning the 528 B cells in the subsampled data to be causal (instead of the 500 randomly selected cells, varying effect size while fixing 25% overlap). The FDR and power are based on applying the Benjamini-Hochberg procedure 65 Table 1). All diseases/traits were well powered (heritability z-score >5), except celiac disease (celiac), systemic lupus erythematosus, MS, subject well-being and type 1 diabetes, which were included due to their clinical importance. The major histocompatibility complex region was removed from all analyses because of its unusual linkage disequilibrium and genetic architecture 67 .
We analyzed 16 scRNA-seq or snRNA-seq data sets. The three atlas-level data sets (TMS FACS, TMS droplet and TS FACS) allow us to broadly associate diverse cell types/populations with disease and compare results between species (mouse/ human) and between technologies (FACS/droplet). The other 13 data sets focus on a single tissue and contain finer-grained annotations of cell types/states and/ or experimentally determined annotations, which allow for better validation (Supplementary Table 2).

Analysis of T cells and autoimmune diseases.
We collectively analyzed all TMS FACS T cells (4,125 cells labeled as "CD4 + α-β T cell", "CD8 + α-β T cell", "regulatory T cell", "mature NK T cell", "mature α-β T cell", or "T cell" in the TMS data; Supplementary Table 5); the more general terms like "T cell" and "mature α-β T cell" were used for cells whose more specific identities were not clear. We processed the T cells following the same procedure as described in the original paper 17,64 . First, we performed size factor normalization (10,000 counts per cell) and log transformation. Second, we selected highly variable genes and computed the batch-corrected principal-component analysis embedding using Harmony 68 , treating each mouse as a batch. Finally, we constructed k-nearest neighbor graphs and clustered the cells using the Leiden algorithm 69 (resolution = 0.7), followed by computing the UMAP embedding. We removed 376 cells from small clusters (less than 100 cells) or whose identities are ambiguous, resulting in 3,769 cells. We annotated the clusters based on the major TMS cell types in the cluster; the label "mature α-β T cell" was omitted because a more specific TMS cell-type label (for example, "CD8 + α-β T") was available in the corresponding cluster. We further characterized disease-associated T cell subpopulations based on marker gene expression, automatic T cell subtype annotation 70 and overlap of specifically expressed genes in each subpopulation with T cell signature gene sets (Supplementary Figs. 14-17 and Supplementary Note). We considered cells from clusters 1-4 as clear CD4 + T cells (1,686 cells) and cells from clusters 1, 2 and 7-9 as clear CD8 + T cells (2,197 cells; the shared clusters 1 and 2 contain a mix of naive CD4 + and CD8 + T cells). We used diffusion pseudotime 38 to assign effectorness gradient for CD4 + and CD8 + T cells separately, where we used the leftmost cell in cluster 2 on the UMAP as the root cell (clearly naive T cell).
We used MSigDB 71 (v7.1) to curate T cell signature gene sets, including naive CD4, memory CD4, effector CD4, naive CD8, memory CD8, effector CD8, T reg , Th1 (T helper 1), Th2 (T helper 2) and Th17 (T helper 17) signatures. For each T cell signature gene set, we identified a set of relevant MSigDB gene sets (22-34  gene sets; Supplementary Table 17), followed by selecting the top 100 most frequent genes in these MSigDB gene sets as the T cell signature genes; a gene was required to appear at least twice and genes appearing the same number of times were all included, resulting in 62 to 513 genes for the 10 T cell signature gene sets (Supplementary Table 24). For gold-standard gene sets used in the analysis of disease gene prioritization, we curated 27 putative drug target gene sets from Open Targets 41 (mapped to 27 of the 74 diseases/traits); for a given disease, we selected all genes with drug score >0 (clinical trial phase 1 and above) and only considered diseases with at least 10 putative drug target genes. We curated 16 Mendelian diseases gene sets from Freund et al. 42 (mapped to 45 of the 74 diseases/traits) (Supplementary Table 21). For comparison of two gene sets, the P value is based on two-sided Fisher's exact test and excess overlap is defined as the ratio between the observed overlap of the two gene sets and the expected overlap (by chance). Of note, for a given query gene set with a fixed size and a fixed level of excess overlap with the reference gene set, the −log 10 P value increases with the size of the reference gene set; we report both excess overlap and −log 10 P value while using the former as our primary metric, which is more interpretable.
Analysis of neurons and brain-related diseases/traits. For the TMS FACS data, we focused on the 484 neurons (TMS label "neuron", excluding cells with TMS label "medium spiny neuron" or "interneuron"). For the Zeisel & Munõz-Manchado et al. data, we applied scDRS to all 3,005 cells and then focused on the 827 CA1 pyramidal neurons ("level1class" label "pyramidal CA1"). For inferring spatial coordinates, we curated differentially expressed genes for each of the six spatial regions (dorsal versus ventral, ventral versus dorsal, proximal versus distal, distal versus proximal, deep versus superficial and superficial versus deep) using the gene expression data from Cembrowski et al. 48 (GEO GSE67403; gene sets in Supplementary Table 24). For each differential gene expression analysis, we selected genes based on fragments per kilobase of transcript per million mapped reads (FPKM) >10 for the average expression in the enriched region (for example, dorsal for the dorsal versus ventral comparison), q-value < 0.05, and log 2 (fold change) >2. We used scDRS and these signature gene sets to assign six spatial scores for each cell. For the regression analysis, we separately regressed the scDRS disease scores for each of the seven brain-related traits (and height, a negative control trait) on each of the six spatial scores. We performed marginal regression instead of joint regression for these spatial scores because the inferred spatial scores for opposite regions on the same axis (for example, dorsal versus ventral) were highly collinear (strongly negatively correlated), and the inferred spatial scores for dorsal, proximal and deep regions (which had strong marginal associations to diseases) had very low pairwise correlations (average |r| = 0.10; Supplementary  Fig. 23d), suggesting these associations were independent.

Analysis of hepatocytes and metabolic traits.
We considered all hepatocytes in the TMS FACS data (1,162 cells). Because the original study clustered all cells from the liver together 17 (limiting the resolution for distinguishing cell states within hepatocytes), we reprocessed and reclustered these cells following the same procedure as we did for the T cells. We further filtered out low-quality cells (proportion of mitochondrial gene counts ≥0.3; likely apoptotic or lysing), resulting in 1,102 hepatocytes. We computed the polyploidy, pericentral and periportal scores by applying scDRS to published polyploidy/zonation signature gene sets (instead of MAGMA putative disease gene sets). We curated signature gene sets for ploidy level, zonation (pericentral/periportal) and putative zonated pathways. We curated four sets of polyploidy signatures, including differentially expressed genes for partial hepatectomy (PH) versus pre-PH 58 (used for the polyploidy score), Cdk1 knockout (case) versus control 58 , 4n versus 2n hepatocytes 60 , large versus small hepatocytes 58 . We curated three sets of diploidy signatures, including differentially expressed genes for pre-PH versus PH 58 , control versus Cdk1 knockout 58 and 2n versus 4n hepatocytes 60 . We curated signature gene sets for pericentral and periportal hepatocytes from Halpern et al. 59 . We curated gene sets for putative zonated pathways from MSigDB 71 (v7.1), including glycolysis (pericentral), bile acid production (pericentral), lipogenesis (pericentral), xenobiotic metabolism (pericentral), beta oxidation (periportal), cholesterol biosynthesis (periportal), protein secretion (periportal) and gluconeogenesis (periportal) (Supplementary Table 24). For the joint regression analysis of scDRS disease score on ploidy and zonation scores, we regressed the polyploidy score out of both the pericentral and periportal score before the joint regression, because the ploidy level confounded both zonation scores. We performed joint regression instead of marginal regression here (unlike the regression analysis in the neuron section), because the polyploidy score was positively correlated with the pericentral and periportal scores (unlike the analysis in the neuron section where the three sets of scores had low correlations).

Statistics and reproducibility.
We analyzed only existing data sets. No statistical method was used to predetermine sample size. No data were excluded from the analyses. We did not use any study design that required randomization or blinding. We replicated our results by performing the same analyses on additional independent data sets; all attempts at replication were successful. Fig. 1 | Additional null simulations. We performed null simulations for various numbers of putative disease genes (100, 500, 1,000, and 2,000 for the four columns respectively) and various types of genes to randomly sample from: all genes (first row), and top 25% genes with high expression (second row), top 25% genes with high expression variance (third row), top 25% overdispersed genes (fourth row). We considered two additional versions of scDRS: scDRS-bin-gs (binary gene sets instead of MAGMA z-score gene weights) and scDRS-adj-ctp (adjusting for cell-type proportion). For scDRS-adj-ctp, we simulated random biased gene sets (high-mean/high-variance/overdispersed) based on the balanced data (inversely weighting cells by cell-type proportion) to better match the model assumption, namely testing for excess expression relative to cells in the balanced data. In each panel, the x-axis denotes theoretical −log 10 p-value quantiles and the y-axis denotes actual −log 10 p-value quantiles for different methods. The 3 versions of scDRS produced well-calibrated p-values in most settings and suffered slightly inflated type I error in panels o,p, possibly because it is hard to match a large number of overdispersed putative disease genes using the remaining set of genes. In comparison, all other methods are less well-calibrated and are particularly problematic when the numbers of putative disease genes are small. Error bars denote 95% confidence intervals around the mean of 100 simulation replicates. Fig. 3 | complete results for cell type-level disease associations for 74 diseases/traits and TMS FAcS 120 cell types. Each row represents a disease/trait and each column represents a cell type (number of cells in parentheses). Heatmap colors denote the proportion of significantly associated cells (FDR < 0.1 across all cells for a given disease). Squares denote significant cell-type-disease associations (FDR < 0.05 across all pairs of the 120 cell types and 74 diseases/traits; 597 significant pairs; MC test; Methods). Cross symbols denote significant heterogeneity in association with disease across individual cells within a given cell type (FDR < 0.05 across all pairs; 273 significant pairs; MC test; Methods). Heatmap colors and cross symbols are omitted for cell-type-disease pairs with nonsignificant cell-type-disease associations. Within the blood/immune block (40 cell types and 21 diseases/traits), 136 of 264 cell-type-disease pairs with significant association also had significant heterogeneity. Within the brain block (11 cell types and 21 diseases/traits), 64 of 133 cell-type-disease pairs with significant association also had significant heterogeneity. Within the other block (69 cell types and 32 diseases/traits), 54 of 146 cell-type-disease pairs with significant association also had significant heterogeneity. We discuss the results for FEV1/FVC. We identified 20 cell types associated with FEV1/FVC (FDR < 0.05), including 5 lung cell types and 15 cell types from other tissues. They can be categorized into 5 sets of associations: (1) type II pneumocyte (2) skin-related cells (3) smooth muscle cells (4) fibroblast-and-MSC-like cells (5) pericyte-like cells. The first 4 sets of associations are consistent with a previous work 75 . The 5th set of pericyte associations is also plausible because pericytes are known to regulate lung morphogenesis 76 . We note that the cell type associations from the lung are more likely to be causal and those from the other tissues are more likely tagging the causal cell types due to shared expression. Numerical results are reported in Supplementary Table 12.

Extended Data Fig. 4 | comparison of cell type-level disease association results between TMS FAcS and TMS droplet (different technologies), TS FAcS (different species).
(a-c) Results for disease association at the cell type-level for TMS FACS, TMS droplet, and TS FACS for diseases and cell types in the blood/immune block (upper left) and the other cell types/diseases block (lower right) in Fig. 3 (TMS droplet and TS FACS do not contain brain data; Supplementary Tables 6,7). The plotting style is the same as Fig. 3. Heatmap colors for each cell-type-disease pair denote the proportion of significantly associated cells (FDR < 0.1); squares denote significant cell-type-disease associations (FDR < 0.05); and cross symbols denote significant heterogeneity in association with disease across individual cells within a given cell type (FDR < 0.05). Heatmap colors (>10% of cells associated) and cross symbols are omitted for cell-type-disease pairs with nonsignificant cell-type-disease associations via MC test. We matched each TMS FACS cell type using the closest cell type in the TMS droplet and TS FACS data; unmatched cell types were colored in grey. (d) Overlap of significant cell-type-disease associations between TMS FACS and TMS droplet (P = 2.8 × 10 −24 , two-sided Fisher's exact test). (e) Pearson's correlation of −log 10 p-values for cell-typedisease associations between TMS FACS and TMS droplet. (f) Overlap of significant cell-type-disease associations between TMS FACS and TS FACS (P = 1.3 × 10 −7 , two-sided Fisher's exact test). (g) Pearson's correlation of −log 10 p-values for cell-type-disease associations between TMS FACS and TS FACS. We determined that the results are highly consistent between TMS FACS and TMS droplet, and are reasonably consistent between TMS FACS and TS FACS. Our method is underpowered in the TS FACS data, possibly due to the smaller sample size (27K cells in TS FACS versus 110K cells in TMS FACS). The current TS FACS data corresponds to the initial data release and there will likely be more cells in future releases 20 .

Extended Data Fig. 5 | Optimizing parameters of scDRS based on expected and unexpected control cell types across 20 traits.
We considered different versions of scDRS by varying methods for selecting (1) putative disease genes (2) weights for the disease genes (3) MAGMA window size. We considered 6 methods for selecting putative disease genes, 4 methods for selecting gene weights, and 3 MAGMA gene window sizes (Supplementary Note). We applied each version of scDRS to the subsampled TMS FACS data (20 repetitions with 10K cells each) and a curated set of 20 traits with expected and unexpected disease-critical cell types (Supplementary Table 15). For a given scDRS version and a given trait, we computed the t-statistic between cells from the expected and unexpected cell types, and divided it by the average t-statistics of results of the given trait from all data sets and all scDRS versions to correct for trait-specific baseline. We evaluated each version by first computing the mean and SE of the normalized t-statistics for a given trait across the 20 repetitions and then combining the estimates across the 20 traits via random-effect meta-analysis. We compared the performance of a pair of scDRS versions by applying the same procedure to the difference of the normalized t-statistics between the two versions. (a) Varying gene selection methods while fixing other parameters as the default. (b) Varying gene weighting methods while fixing other parameters as the default. (c) Varying MAGMA gene window size while fixing other parameters as the default. The default version was denoted in dark blue. Error bars denote 95% confidence intervals around the mean based on meta-analysis across 20 subsampled data sets and 20 traits, using procedures as described above. * denotes P < 0.05 and ** denotes P < 0.001 for significant differences relative to the default configuration; one-sided tests based on the estimated mean and CIs. Numerical results are reported in Supplementary Table 16. Fig. 8 | complete results of correlations between scDRS disease scores and inferred spatial coordinates across cA1 pyramidal neurons in 7 single-cell data sets (extending results in Fig. 5b). (a) Results for regressing the scDRS disease scores against the inferred spatial coordinates for each disease/trait and each inferred spatial coordinate. Color represents the t-statistics and stars represent significant associations (* denotes P < 0.05 and ** denotes P < 0.005, one-sided MC test; Methods). Heatmap color represent the average t-statistics across the 7 brain-related diseases/traits (excluding height) for each data set and stars represent significant associations by combining p-values across data sets using Fisher's combined probability test. (c) Summary of the association between brain-related diseases and the inferred spatial coordinates for the mouse and human data sets in panel b.