The open GWAS database¹ was utilized to obtain GWAS statistics for different types of stroke, including ebi-a-GCST006906 for stroke, ebi-a-GCST005840 for large artery atherosclerosis stroke, and ebi-a-GCST006909 for small vessel stroke (Malik et al., 2018). To ensure adequate overlap between the genes obtained from single-cell data and the GWAS risk loci, we integrated a brain vascular single-cell atlas (Vanlandewijck et al., 2018), which was obtained using SMART-seq2-based technology. Specifically, we analyzed the brain vascular single-cell dataset from the Gene Expression Omnibus database (GEO; https://www.ncbi.nlm.nih.gov/geo/) using accession numbers GSE98816 (Vanlandewijck et al., 2018).

We processed and integrated our raw single-cell count matrices using the Seurat package (version 4.0; https://satijalab.org/seurat/). To ensure data quality, we excluded low-quality cells which expressed ≤200 genes, ≥20,000 unique molecular identifiers, or ≥20% mitochondrial genes. After that, we performed normalization and dimensionality reduction using the t-Distributed Stochastic Neighbor Embedding (t-SNE) algorithm (Skrodzki et al., 2024) for non-linear dimensionality reduction. Clusters were identified and annotated using previously published data (He et al., 2018).

To identify a set of genes that are expressed in specific cell types and to quantify their expression specificity (ES), we used CELLEX (CELL type EXpression-specificity) (Timshel et al., 2020). This tool allowed us to measure the ES of each gene in the expression matrix of each cell type. To simplify the analysis, we calculated a single ES estimate (ESμ) for each gene based on its expression in each cell type. This score represents the degree of specificity of the gene’s expression to the cell type in which it is expressed.

To systematically prioritize cell types that may play etiological roles in complex diseases such as stroke, we employed an integrative analytical framework that combines scRNA-seq data with GWAS summary statistics. Specifically, we utilized the CELLECT tool (Timshel et al., 2020), which incorporates gene expression specificity information across cell types to evaluate heritability enrichment using S-LDSC. The CELLECT workflow consists of the following key steps: (1) Calculation of cell type-specific expression profiles: For each cell type or cluster identified in the scRNA-seq dataset, CELLECT quantifies expression specificity by computing metrics such as the mean expression and ESμ score (expression specificity metric), identifying genes that are highly specific to each cell Type. (2) Gene annotation to regulatory regions: Each cell type-specific gene is annotated to a genomic region that includes a ± 100 kilobase (kb) window surrounding the transcription start site, capturing nearby regulatory elements that may influence gene expression. (3) Linkage disequilibrium (LD) score computation: For each annotated gene region, CELLECT calculates LD scores for HapMap3 SNPs, quantifying the extent to which these variants are correlated with each other across the population. This step creates a cell type-specific annotation map that reflects the distribution of regulatory elements across the genome. (4) Heritability enrichment analysis: Finally, the framework assesses whether the SNPs associated with a given trait or disease (e.g., stroke or stroke subtypes) are significantly enriched in the gene sets specific to certain cell types. This is done through S-LDSC, which models the contribution of each annotated SNP set to the overall SNP-based heritability of the trait. A statistically significant enrichment indicates that the corresponding cell type may be critically involved in the pathogenesis of the disease. By integrating GWAS summary statistics for stroke and its subtypes with high-resolution scRNA-seq data from cerebrovascular cells, this approach enables the identification of cell types whose gene expression landscapes are disproportionately enriched for genetic risk variants, thereby pinpointing potential cellular origins of disease susceptibility.

To identify differentially expressed genes between cell subpopulations, we utilized the “Findmarker” function in the Seurat package. GSEA was performed using the “clusterProfiler” package.

We analyzed the developmental trajectory of cell types using the Monocle2 (v2.12.0) R package (Trapnell et al., 2014; Qiu et al., 2017). To identify gene ordering (qval <0.01), we utilized the differential GeneTest function of the Monocle2 package. Subsequently, we utilized the reduce Dimension function and performed the order Cells function with default parameters to infer the trajectory.

pySCIENIC (Single-Cell rEgulatory Network Inference and Clustering)(version 0.11.0) was used to identify single-cell gene regulatory network (GRN) (Van de Sande et al., 2020). To construct the GRNs, we first generated an expression matrix and co-expression modules of transcription factors (TFs) and their potential target genes based on co-expression analysis. We then performed TF motif analysis using RcisTarget to identify direct-binding regulons. To evaluate the activity of the regulons in each cell, we used AUCell (regulon_AUC) to calculate the regulon activity score (RAS), which represents the average activity score across regulons. Finally, we summarized the GRNs and regulons for the endothelial cell types of mouse cerebrovascular cells.

We conducted a comprehensive analysis of a scRNA-seq dataset obtained from a recent study (He et al., 2018), which aimed to uncover the gene expression profiles of vascular and vessel-associated cells across defined cerebrovascular zonation in the mouse brain. Our analysis identified 14 distinct cell clusters using t-Distributed stochastic neighbor embedding (Figure 1A), and revealed the top genes that were specifically expressed in each cell cluster (Figure 1B). Notably, genes such as Cldn5, Myh11, Pla1a, Theme119, Plp1, and Lum, were consistently expressed at higher levels in their respective cell clusters (Figure 1C).

To quantitatively assess and rank gene expression specificity in the 14 cerebrovascular cell subtypes, we utilized the CELLEX toolkit. The ESμ score was utilized to represent the degree of gene expression specificity in each cell cluster. Our analysis identified the most specifically expressed genes in each cell cluster (Figure 1D), and their expression patterns were consistent with those observed in the dot plot of average expression. These findings indirectly validated the reliability of the ESμ score in assessing gene expression specificity.

To identify vascular cell types associated with stroke, we employed the CELLECT algorithm to analyze the 14 vascular cell types using three different GWAS statistics, namely stroke, Stroke_Large_AS, and ischemic Stroke_small_vessel. Our analysis revealed that two of the EC clusters, namely the aEC cluster and EC2 cluster, showed significant enrichment for variants associated with Stroke_large_AS (Figure 2A). Moreover, we assessed the ESμof seven stroke candidate genes obtained from a GWAS study (Bevan et al., 2012), which supported the reliability of our findings (Figure 2B). Notably, several of these candidate genes, including DDAH1, NOS3, and SGK1, are specifically enriched in either the aEC_cluster or the EC2_cluster, providing insights into the molecular mechanisms by which these cell types contribute to the pathogenesis of large artery atherosclerotic stroke.

To explore the relationship between stroke pathology and the two distinct subtypes of endothelial cells, we employed GSEA to analyze their cellular functions (Figures 3A–D). Our findings revealed significant enrichment of specific pathways in each subtype. In the aEC_cluster, we identified a significant upregulation of pathways associated with fluid shear stress and atherosclerosis, growth hormone synthesis, inflammatory mediator regulation of TRP channels, and TNF signaling (Figure 3C). Conversely, the EC2_cluster exhibited a significant downregulation of pathways related to muscle organ development, muscle tissue development, and regulation of epithelial cell proliferation (Figure 3D). These findings suggest distinct functional profiles and underlying molecular processes associated with the two critical endothelial cell subtypes.

Moreover, we identified several key molecules, including Nr4a1, Fos, and Alox12, that exhibited significant changes in both aEC_cluster and EC2_cluster (Figure 3E). These molecules are known to be involved in crucial biological processes, highlighting their potential roles in mediating the pathogenesis of stroke within these endothelial cell subtypes.

To gain insights into the differentiation trajectory of the two stroke-associated endothelial cell subtypes, we performed pseudotime analysis on all endothelial cells using Monocle2. Our analysis revealed a differentiation trajectory characterized by the bifurcation of endothelial cells into two states: state 2 cells differentiating toward either state 1 or state 3. Specifically, the aEC_Cluster mainly belonged to state 1, while EC2_Cluster predominantly belonged to state 2. As a result, the two pathogenic EC classes within the EC2_Cluster gradually transitioned toward aEC_Cluster differentiation (Figures 4A–C). Notably, the aEC_Cluster represents the terminal stage of differentiation.

Furthermore, we performed unsupervised clustering of gene expression kinetic patterns and identified six expression modules shared by endothelial cells. Intriguingly, genes within modules 2, 3 and 5 showed a gradual elevation in expression as pseudotime advanced (Figure 4D). Consistent with our findings in Figure 3, differential expression and functional enrichment analysis revealed the gradual upregulation of genes associated with atherosclerosis, such as FOS, NR4A1, and alox12, along the differentiation trajectory from EC2_Cluster to aEC_Cluster.

To investigate the upstream transcriptional regulatory events that play potential key roles in the differentiation of endothelial cells into these critical cell subtypes, we conducted a gene regulatory network analysis. Using SCENIC’s integrated algorithm, we identified 10 transcriptional regulator modules of regulons (i.e., transcription factor with its target) and their co-expression modules in all the endothelial cells (Figure 5A; Supplementary Table 1). We further evaluated the scoring of each endothelial cell subtype in these modules and found that aEC_cluster had the highest scores in modules M1 and M2, while EC2_cluster had the highest score in M4 (Supplementary Figure 1). To identify the key transcriptional regulatory network in the two potentially critical cell subtypes associated with stroke, we first scored the regulons in each subtype (Supplementary Table 2). From these results, the top regulons enriched in each subtype were identified (Figure 5B). Notably, the results revealed that the activity of the top five regulons, namely TF_Mafk, TF_Gata2, TF_Irf7, TF_Stat6 and TF_Jund, was significantly elevated in aEC (Figure 5B). Further analysis revealed the extensive regulation of molecular networks in aEC by transcription factors Gata2, Irf7, and Jund. Remarkably, this regulation included molecules such as Fos and Nr4a1, which were displayed progressively elevation during aEC_cluster differentiation in the pseudotime analysis (Figures 4D, 5C). To further validate the reliability of our transcriptional regulatory factor predictions, we incorporated scRNA-seq data from an external model of middle cerebral artery occlusion (MCAO) in the brain (Zhang et al., 2023). The results indicated that the levels of Ifr7 and Fos in the MCAO model group were significantly higher than those in the control group (Figure 5D), thereby suggesting the reliability of our findings to some extent.