hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data

1The development and widespread adoption of single-cell and spatial genomics approaches has led to routine generation of high-dimensional datasets in a variety of biological systems. These technologies are frequently used to study developmental stages, evolutionary trajectories, disease states, drug perturbations, and other experimental conditions. Despite the inherent complexity and interconnectedness of biological systems, studies leveraging single-cell and spatial genomics typically analyze individual features (genes, isoforms, proteins, etc.) one by one, greatly oversimplifying the underlying biology. These datasets provide an opportunity for investigating and quantifying the relationships between these features to further contextualize their roles across biological conditions of interest.

2Here we developed hdWGCNA, a framework for co-expression network analysis in single-cell and spatial transcriptomics data. Co-expression networks are based on transformed pairwise correlations of input features, resulting in a quantitative measure of relatedness between genes., Hierarchical clustering on the network structure allows us to uncover functional modules of genes whose expression profiles are tightly intertwined,, which typically correspond to specific biological processes and disease states. Considering that unique cell types and cell states have distinct gene expression programs, we designed hdWGCNA to facilitate multi-scale analysis of cellular and spatial hierarchies. hdWGCNA provides a rich suite of functions for data analysis and visualization, providing biological context for co-expression networks by leveraging a variety of biological knowledge databases. To maximize usability among the genomics community, the hdWGCNA R package extends the data structures and functionality of the widely used Seurat package,,, and we developed an extensive documentation website for hdWGCNA demonstrating its use on new datasets. Further, we used hdWGCNA to analyze a single-cell RNA sequencing (scRNA-seq) dataset consisting of 1 million cells, showcasing the scalability of hdWGCNA in large datasets.

3In this study, we applied hdWGCNA in a variety of high-dimensional transcriptomics datasets from different technologies and biological conditions. As a common use case, we first performed iterative network analysis of the major cell types in the human prefrontal cortex (PFC), identifying shared and specific network modules in each cell type. We constructed co-expression networks in anterior and posterior mouse brain sections profiled with Genomics Visium spatial transcriptomics (ST), and found distinct spatial patterns of these gene expression programs. Using long-read scRNA-seq data from the mouse hippocampus, we uncovered splicing isoform co-expression networks in the radial glia lineage involved in cell fate specification. Network analysis of inhibitory neurons from published single-nucleus RNA sequencing (snRNA-seq) in autism spectrum disorder (ASD) donors revealed modules disrupted in ASD containing key genetic risk genes such as SCN2A, TSC1, and SHANK2. We performed consensus co-expression network analysis of microglia from three Alzheimer’s disease (AD) snRNA-seq studies,,, yielding multiple gene modules corresponding to disease-associated microglia and polygenic risk of AD. Finally, we used hdWGCNA to project gene modules from two bulk RNA-seq studies of AD patients into an snRNA-seq dataset of the AD brain, showing that our approach allows for interrogation of gene modules and networks that have been previously identified.

4Here we describe hdWGCNA, a comprehensive framework for constructing and analyzing co-expression networks in high-dimensional transcriptomic data (Figure 1A). Given a gene expression dataset as input, co-expression network analysis typically consists of the following analysis steps: computing pairwise correlations of input features, weighting correlations with a soft-power threshold, computing the topological overlap between features, and unsupervised clustering via the Dynamic Tree Cut algorithm (Figure S1 and STAR Methods). The sparsity and noise inherent in single-cell data can lead to spurious gene-gene correlations, thereby complicating co-expression network analysis. Additionally, the correlation structure of single-cell or spatial transcriptomic data varies greatly for different subsets (cell types, cell states, anatomical regions). A typical hdWGCNA workflow in scRNA-seq data accounts for these considerations by collapsing highly similar cells into “metacells” to reduce sparsity while retaining cellular heterogeneity and by allowing for a modular design to perform separate network analyses in specified cell populations.

5Metacells are defined as small groups of transcriptomically similar cells representing distinctive cell states. There are several approaches to identify metacells from single-cell genomics data.,,, We leverage a bootstrapped aggregation (bagging) algorithm for constructing metacell transcriptomic profiles from single-cell datasets by applying K-nearest neighbors (KNN) to a dimensionality-reduced representation of the input dataset (STAR Methods, Algorithm 1). This approach can be performed for each biological replicate to ensure that critical information about each sample (age, sex, disease status, etc.) is retained for downstream analysis. We computed gene-gene correlations in the normalized gene expression matrix from the single-cell dataset and metacell expression matrices while varying the number of cells to collapse into a single metacell (the KNN parameter). The distribution of these gene-gene correlations displays a spike at zero for the single-cell expression matrix, with flattened distributions corresponding to more non-zero correlations in the metacell matrices, indicating that metacell expression profiles are less prone to noisy gene-gene correlations compared with the single-cell matrix (Figure 1B) (STAR Methods). We note that sparsity (defined in Equation 1) is greatly reduced in the metacell matrices for each cell type compared with the single-cell matrices, with over a 10-fold reduction in some cases (Figure 1C). We applied hdWGCNA to a dataset of CD34+ hematopoietic stem and progenitor stem cells using two additional metacell approaches, and found that all approaches were suitable for downstream network analysis (Figure S2; Data S1). Metacell algorithms strive to retain biologically meaningful signals spanning a spectrum of cell states in a tissue of interest; therefore, it is necessary to carefully apply these approaches to avoid obscuring these cell states. For example, the hdWGCNA metacell algorithm requires a dimensional reduction of the input expression matrix, but these reductions often contain technical artifacts. The choice of dimensionality reduction method and handling of technical artifacts would then influence the effectiveness of metacell construction. Further, the optimal number of cells to merge together to form a single metacell may differ across cell types and tissues, attempting to balance between increasing information content of the aggregated group while avoiding merging of dissimilar cells. Aside from metacell approaches, pseudo-bulk aggregation of all cells in a given population have yielded favorable results in benchmarks of differential gene expression tests, suggesting that, given a sufficient sample size, pseudo-bulk expression profiles are likely suitable for co-expression network analysis.

6While co-expression modules consist of many genes, it is convenient to summarize the expression of the entire module into a single metric. Module eigengenes (MEs), defined as the first principal component of the module’s gene expression matrix (STAR Methods, Algorithm 2), describe the expression patterns of entire co-expression modules. hdWGCNA computes MEs using specific accommodations for high-dimensional data, allowing for batch correction and regression of continuous covariates (STAR Methods, Algorithm 2). Optionally, hdWGCNA can use alternative gene scoring methods such as or UCell or Seurat’s AddModuleScore function, and we show that these scores are correlated with MEs (Figure S3).

7We demonstrate hdWGCNA in single-cell transcriptomic data through an iterative network analysis of six major cell types in the Zhou et al. human PFC snRNA-seq dataset of 11 cognitively normal donors (Figure 1A). We constructed metacells and performed co-expression network analysis for each major cell type in the human PFC dataset using the standard hdWGCNA workflow, yielding distinct network structures and sets of gene modules (Data S2 and S3). Networks were constructed using metacell expression matrices for each cell type separately, but we computed MEs for each module using the entire snRNA-seq dataset, allowing us to interrogate the cell-type specificity of these modules’ expression programs across all cell types. This iterative network analysis revealed 96 co-expression modules across the six major cell types. Through differential module eigengene (DME) analysis, we found shared and distinct module expression patterns across different cell types (Data S2; STAR Methods), and we highlight specific modules from each cell type (Figures 1D and 1E). Further, we performed a pairwise gene set overlap analysis of the 96 co-expression modules, and, while we did find that some modules had significant overlaps across the different cell types, the gene sets comprising these modules were overall quite distinct, with a maximum Jaccard index between two modules of 0.297 and a median of 0.005 (STAR Methods and Figure S4). The expression of module hub genes, which are highly connected members of the co-expression network ranked by eigengene-based connectivity (kME), tend to display cell-type-specific patterns, such as the myelination genes CNP and PLP1 in oligodendrocyte (ODC) module ODC-M3 (Figure 1D). However, some co-expression modules may correspond to cellular processes common to multiple cell types, in which case the hub genes may be widely expressed. We inspected the MEs of selected cell-type-specific modules and found that the overall expression patterns were similar to that of their constituent hub genes (Figures 1D and 1E).

8We showcase some of the downstream functionalities of hdWGCNA using the ODC co-expression network (Figures 1F–1H). For network visualization, we used Uniform Manifold Approximation and Projection (UMAP) to embed the co-expression network topological overlap matrix (TOM) into a two-dimensional manifold, using the topological overlap of each gene with the top hub genes from each module as input features (STAR Methods; Figure 1F). We found that eight of the 10 ODC modules were specifically expressed in ODC cells based on their MEs (Figure 1G; Wilcoxon rank-sum test Bonferroni-adjusted p <0.05). Finally, we performed module preservation analysis to test the reproducibility of these modules in an independent dataset and found that all of the ODC-specific modules were significantly preserved ( summary preservation ). In sum, these network analyses in the human PFC dataset shows the core capabilities of the hdWGCNA workflow (Figure S1). Finally, we performed a similar iterative network analysis on a peripheral blood mononuclear cell (PBMC) scRNA-seq dataset of nearly 1M cells, highlighting the scalability of hdWGCNA in large datasets (Figure S5; Data S1).

9We measured the runtime and memory usage of hdWGCNA as a function of the number of input cells. Using the 65,415 neuronal cells from the Velmeshev et al. human PFC snRNA-seq dataset (54 samples), we ran hdWGCNA on different-sized subsets ranging from 1,000 to 50,000 cells to test the runtime and memory consumption of the main network analysis steps (Figures 2A and 2B). We report the memory upper bound in gigabytes measured throughout the duration of each function. The runtime of the MetacellsByGroups function increased steadily with the number of cells, but the memory usage plateaus. This function attempts to construct a target number of metacells within each biological replicate and each cell population, and the algorithm terminates early if this target is reached, thus explaining the plateau in the memory usage graph. While TestSoftPowers generally had a low memory footprint, it was the slowest individual function based on these tests. Importantly, TestSoftPowers can be sped up by using a subset of the data, or by testing fewer soft-power thresholds than the default. The efficiency of ConstructNetwork varies both with the number of input cells and features, where this calculation will slow down as more cells and features are included. ModuleEigengenes uses the implicitly restarted Lanczos bidiagonalization algorithm (IRLBA) for fast singular value decomposition (SVD) of sparse matrices, and the runtime and memory usage of this function both linearly increase with the number of cells in the dataset. Optionally, ModuleEigengenes can employ the harmony algorithm following SVD, which increases runtime but not memory usage. Further, the efficiency of this function varies with the number of co-expression modules detected and with the number of features in each modules. Finally, ModuleConnectivity computes eigengene-based connectivity as product-moment correlation coefficients between the sparse gene expression matrix and the MEs matrix, which resulted in fast calculations with low memory usage.

10We next sought to evaluate the co-expression networks identified by hdWGCNA using a functional coherence analysis. We used the EGAD neighbor-voting algorithm to predict known biological pathway associations of genes based on the co-expression network structure using the cell-type-specific co-expression networks from the Zhou et al. human PFC dataset. In principle, we expect that co-expressed genes are involved in similar biological processes, and therefore co-expression network structures should be predictive of biological pathway membership. Briefly, EGAD performs a 3-fold cross-validation classification by occluding the pathway labels of a subset of genes and then attempting to predict the pathway membership of those occluded genes based on their labeled neighbors in the network. We used EGAD to test the functional coherence of our hdWGCNA co-expression networks for a set of Gene Ontology (GO) terms, reporting area under the receiver operating characteristic curve (AUC) value for each term (Figure 2C). We report a similar level of functional coherence in these co-expression networks to a previous study that evaluated co-expression networks derived from scRNA-seq data with different measures of gene-gene association. The inhibitory neuron network performed the best for functional coherence with a median neighbor-voting AUC of 0.592, while the lowest-performing network was from the oligodendrocytes with a median AUC of 0.549. We tested whether there was a bias toward genes that were multifunctional based on the frequency that they appeared in the annotated set of GO terms, and we found that multifunctional genes did not bias the co-expression functional coherence results (Figure 2D).

11In principle, genes within the same co-expression modules derived from specific cell types should be functionally related or co-regulated. The expression of module hub genes, which exhibit the highest intramodular connectivity, may be predictive of the expression of other module member genes if the network is well defined and contains meaningful structures. For each of the 96 cell-type PFC co-expression modules, we sought to predict the expression of each gene using the top 10 module hub genes as the input features to a XGBoost regularized regression model. In this analysis, we performed 5-fold cross-validation, and we report the performance as root-mean-square error (RMSE) of the test set averaged over each fold (Figure 2E). Overall, we found that module hub gene expression was generally predictive of module member gene expression across all modules in the six cell-type co-expression networks, where the module with the best performance had an average test set RMSE of 0.0159 and the module with the worst performance had an average test set RMSE of 0.209 (Figure 2E). This analysis and our functional coherence analysis provide support that hdWGCNA co-expression networks and gene modules capture biologically relevant information in specific cell types.

12ST enables the investigation of biological patterns that might otherwise be hidden in other -omics technologies, such as scRNA-seq or bulk RNA-seq., We used hdWGCNA to identify spatial co-expression network modules in the murine brain using a publicly available Visium transcriptomics dataset from Genomics (Figure 3A). This ST dataset consists of one posterior and one anterior slice originating from a sagittal brain section from a single male mouse at 8 weeks of age. Sequencing-based ST approaches such as Visium yield transcriptome-wide gene expression profiles localized to individual “spots” where a single spot likely contains multiple cells, and this dataset is composed of 2,696 spots in the anterior slice and 3,353 spots in the posterior slice. Data sparsity is also inherent to the current generation of these technologies, therefore we propose a metaspot aggregation approach prior to network analysis (Figure S6). Evenly spaced spots throughout the input ST slide are used as principal spots, with at least one other spot in between two principal spots. The transcriptomes of the principal spots and their direct neighbors are aggregated into metaspot expression profiles, containing at most seven ST spots (Figure S6A). Similar to metacells in scRNA-seq, the sparsity of the metaspot expression matrix was reduced compared with the original ST matrix (Figure S6B), and the distribution of gene-gene correlations in the metaspot expression matrix was less concentrated at zero (Figure S6C). hdWGCNA is capable of processing any number of ST samples in the same co-expression network analysis by constructing metaspots separately for each sample.

13We applied hdWGCNA in the mouse brain Visium dataset, identifying 12 spatial modules (SM1-12; Figure S7; Data S1), and we embedded the co-expression network in two dimensions using UMAP (Figure 3B). DME analysis showed that spatial co-expression modules displayed distinct regional expression profiles based on their MEs (Figure 3C; Data S3), encompassing a wide array of cellular processes such as the myelination module SM1 in the white matter tracts, and synaptic transmission modules SM7, SM9, SM11, and SM12 (Figure S7C; Data S2). For example, DME analysis showed that expression of SM4 was localized to the ventricles and cortical layer 1 near the blood-brain barrier (Figure 3C). Further, the hub genes of SM4 include hemoglobin subunits (Hba-a1, Hba-a2, Hbb-bt), and we show that SM4 was enriched for biological processes associated with brain vasculature (Figures 3B and S7C). We compared these gene modules with cluster marker genes from a whole-mouse-brain snRNA-seq dataset and found significant correspondences, such as the striatum module SM7 and medium spiny neurons (Fisher’s exact test false discovery rate [FDR] <0.05; Figure S7D). Additionally, we performed network analysis on a subset of this dataset containing cortical layers 2–6 (Figure S8), identifying additional fine-grained spatial co-expression modules localized to specific cortical layers (Data S1 and S2).

14Different isoforms of the same gene are often involved in distinct biological processes. Conventional single-cell transcriptomics assays capture information at the gene level, thereby missing much of the biological diversity and regulatory mechanisms that occurs at the isoform level. Emerging long-read sequencing approaches enable us to profile cellular transcriptomes at isoform resolution,,,, thus providing new opportunities to model the relationships between isoforms using co-expression network analysis.

15We used hdWGCNA to perform isoform co-expression network analysis in radial glia lineage cells from the mouse hippocampus at postnatal day 7 (P7) profiled with single-cell isoform RNA sequencing (ScISOrSeq) (Figure 4A; STAR Methods). This dataset contains isoform-level and gene-level expression data from 6,832 nuclei derived from a single mouse hippocampus sample. Radial glia, which share transcriptomic similarities with mature astrocytes, are progenitor cells that give rise to numerous distinct cell fates, including neuronal cells, astrocytes, oligodendrocytes, and ependymal cells., To model this developmental process, we applied Monocle3 pseudotime to 2,190 radial glia lineage cells (Figure 4B). We identified three trajectories corresponding to distinct cell fates, termed the ependymal (EPD) trajectory, astrocyte (ASC) trajectory, and the neural intermediate progenitor cell (NPC) trajectory.

16Isoform co-expression network analysis revealed 11 modules in the radial glia lineage (Figure 4C; Data S1). Of the genes retained for network analysis, 61.5% had a single isoform, 18.2% had multiple isoforms that were all assigned to the same module, and 20.4% had multiple isoforms spread across several modules (Figure 4D). Thus, these network modules capture information about the roles of different isoforms of the same gene in distinct biological processes. We inspected module eigenisoform (MEiso) patterns throughout the developmental lineage, thereby uncovering isoform modules critical for cell fate decisions (Figure 4E; Data S2 and S3). Increased expression of modules RGL-M1 and RGL-M2, which were enriched cilium assembly genes (Figure 4F), was associated with the transition from a radial glia to an ependymal cell state. A steady expression level of module RGL-M5 (glial development, astrocyte differentiation) was found in the transition from radial glia to astrocytes, while a decreased expression of RGL-M5 led to alternative fates. Four modules (RGL-M3, RGL-M8, RGL-M9, and RGL-M11) displayed an increase in expression in the neuronal trajectory, containing genes associated with cellular processes such as non-canonical Wnt signaling, neuronal synaptic plasticity, and RNA splicing (Figure 4F).

17We inspected the isoforms of three selected genes that had hub isoforms in different co-expression modules: Gfap, H3f3b, and Cd9 (Figures 4G–4I). Gfap encodes a key intermediate filament protein in astrocytes that is involved in astrocytic reactivity during central nervous system (CNS) injuries or neurodegeneration, and we found that modules RGL-M4 and RGL-M6 contained hub isoforms of Gfap featuring alternative splicing, alternative transcription start site (TSS) usage, and alternative transcription end site (TES) usage (Figure 4G). Different isoforms of the histone H3.3 subunit gene H3f3b were hubs for modules RGL-M1 and RGL-M11, which were associated with ependymal and neuronal cell fates respectively, suggesting that alternative TES usage in H3f3b plays a role in regulating epigenetic factors in murine hippocampal development (Figure 4F). Cd9 encodes a transmembrane protein and is a known glioblastoma biomarker, and we found subtle differences in the TSS between hub isoforms in modules RGL-M6 and RGL-M9 that we show as a splicing summary graph (Figure 4I), supporting functional changes mediated by small isoform differences.

18Co-expression networks can be interrogated to further understand the molecular phenotypes of complex polygenic diseases in primary human tissue samples. We applied hdWGCNA to 20,249 inhibitory neurons (INHs) from an snRNA-seq dataset of the human PFC in 22 ASD patients, 24 age-matched controls, and eight epilepsy patients (Figures 5A and S10; Data S1). The INH network contained 14 modules, and we show hub genes that have a known association with ASD in the SFARI database on the co-expression UMAP (Figure 5B). The MEs showed that some modules were primarily confined to a single INH cluster (INH-M3, INH-M1) while others were spread across multiple neuronal groups (Figure 5C). Furthermore, DME analysis revealed significant differences between MEs in ASD and control samples for all modules except INH-M4 in at least one INH subpopulation (Figure 5D; Data S3; Wilcoxon rank-sum test Bonferroni-adjusted p <0.05). However, by focusing on the DME results with an absolute average (fold change), we note that many of the largest differences were found in the SST+ inhibitory neuron clusters. Furthermore, three co-expression modules (INH-M11, INH-M13, and INH-M3) were significantly enriched in ASD-associated genes from the SFARI database and the latest genome-wide association study (GWAS) of ASD (Figure 5E), but we note that all of these modules contained several ASD-associated SFARI genes.

19INH-M11 was enriched for genes associated with synaptic transmission, ion transport, glutamate receptor signaling, and nervous system development (Figure 5F; Data S2), and this module was downregulated in ASD for five of the six INH subtypes (Figure 5D). Similarly, INH-M13 was associated with RNA processing (Figure 5F) and was downregulated in ASD in all INH subtypes except PVALB+ neurons (Figure 5D). One of the INH-M13 hub genes is CHD2, whose de novo variants have been identified in individuals with ASD., CHD2 is part of the CHD family of chromatin-modifying proteins and can alter gene expression by modification of chromatin structure. Similarly, rare loss-of-function mutations have been reported in the SCN2A gene, a hub gene of the INH-M11 module. We also find enrichment of several ASD-associated genes such as TSC1 (INH-M8), SMARCA4 (INH-M8), SHANK2 (INH-M4), and CPEB4 (INH-M1), highlighting that these modules are functional and provide new insights into the role of inhibitory neurons in ASD. Finally, we tested for the preservation of these modules in 19,425 inhibitory neurons from an snRNA-seq dataset of the PFC from donors with major depressive disorder (MDD) and controls (34 samples), and we found substantial evidence of preservation across all modules except INH-M1 (Figures S10C–S10E).

20Microglia, the resident immune cells of the brain, are implicated in the pathology and genetic risk of several CNS diseases, including AD.,,, Transcriptomic and epigenomic studies in human tissue and AD mouse models have identified multiple cell states of microglia, representing a spectrum between homeostatic and disease-associated microglia (DAMs).,, Our previous study defined a set of transcription factors, genes, and cis-regulatory elements involved in the shift between homeostatic and DAM cell states in human AD, identifying shared and distinct signatures compared with the DAM signature from 5xFAD mice. Here we sought to expand on previous work by providing a systems-level analysis of gene expression throughout the spectrum of microglia cell states.

21We modeled the cell-state continuum between homeostatic and DAM-like microglia by employing a pseudotime analysis of microglia from three human AD snRNA-seq datasets,, (Figures 6A, 6B, and S11). Next, we performed consensus co-expression network analysis using microglia integrated from three human AD snRNA-seq datasets,,, identifying four consensus modules (Figure 6C; Data S1). Consensus network analysis is an approach that performs network analysis separately for each dataset, followed by a procedure to retain structures common across the individual networks, and thus it is well suited for analyzing microglia co-expression from these different sources (STAR Methods).

22Classical markers of homeostatic microglia, such as CSF1R, CX3CR1, and P2RY12, were members of MG-M2, while known DAM genes, including APOE, TYROBP, and B2M, were members of MG-M1. GO term enrichment analysis associated MG-M2 with homeostatic microglia functions such as cell migration, synapse organization, and response to colony-stimulating factor, contrasting disease-related processes enriched in MG-M1 including amyloid fibril formation, microglial activation, maintenance of blood-brain barrier, and cytokine production (Figure 6D; Data S2). Together, this suggests that MG-M1 comprises the gene network underlying DAM activation in AD, while MG-M2 represents the network of homeostatic microglia genes. The MEs for MG-M1 and MG-M2 display opposing patterns throughout the microglia pseudotime trajectory, contextualizing this trajectory as the transcriptional shift from homeostatic microglia (start) to a DAM-like cell state (end) (Figures 6E and 6F). Furthermore, DME analysis revealed significant changes in these modules between AD and control brains in evenly spaced windows throughout the microglia trajectory (Figure 6G; Wilcoxon rank-sum test Bonferroni-adjusted p <0.05; Data S3). Co-expression networks behave as functional biological units; therefore, we reason that the hub genes and other members of MG-M1 represent candidates for an expanded set of human DAM genes including ACTB, TPT1, and EEF1A1.

23Aside from modules MG-M1 and MG-M2, which contained well-known microglia gene signatures, we also identified modules MG-M3 and MG-M4 containing genes associated with key microglial processes such as axon guidance, phagocytosis, and myeloid cell differentiation (Figures 6C and 6D). CD163, a hub gene of MG-M4, is known to be involved in the breakdown of the blood-brain barrier., The trajectory of MG-M4, containing CD163 as a hub gene, was consistent with that of DAM-like module MG-M1, and was enriched for processes including phagocytosis, myeloid cell differentiation, and neutrophil activation (Figure 6D); therefore, it is possible that MG-M4 represents an alternative microglial activation module. We performed single-cell polygenic risk enrichment for AD risk in the microglia trajectory,, and identified a significant increase throughout the trajectory, revealing an enrichment of AD genetic risk single-nucleotide polymorphisms (SNPs) in DAMs (Figure 6H; STAR Methods; Data S3). We show that expression of these modules was significantly correlated with AD genetic risk (Pearson correlation p <0.05), with the strongest correlation in alternative activation module MG-M4 (Figure 6I).

24To ensure that these microglial modules were reproducible across other datasets and in mouse models of AD, we performed module preservation analysis (STAR Methods; Figure 6J). We projected the microglial consensus modules into a dataset of the PFC in aged human samples, the superior frontal gyrus (SFG) and entorhinal cortex (EC) in AD samples, the occipital cortex (OC) and the occipitotemporal cortex (OTC) in human AD samples, the PFC from 5xFAD mice, and whole-brain samples from 5xFAD mice (Figure 6H). Additionally, we projected these modules into an snATAC-seq dataset of the PFC in human AD, using gene activity as a proxy for gene expression from chromatin accessibility data. These module preservation tests showed the microglia consensus modules were broadly preserved and reproducible across brain regions and in mouse models of AD, providing further support that this network is relevant in AD biology and microglial activation.

25hdWGCNA allows for interrogating co-expression modules inferred from a given reference dataset in a query dataset. Modules can be projected across datasets by computing MEs in the query dataset, and preservation of the network structure can be assessed via statistical testing. For example, modules can be projected between different species to link transcriptomic changes between mouse models and human disease patients, or modules can be projected across data modalities from single-cell to spatial transcriptomics to provide regional context to cellular niches.

26To date, it remains cost-prohibitive for most researchers to perform high-dimensional -omics studies of large patient cohorts, but there are numerous large-scale disease-relevant bulk RNA-seq datasets containing thousands of samples from consortia such as the Encyclopedia of DNA Elements (ENCODE), the Genotype-Tissue Expression (GTEx) project, and The Cancer Genome Atlas (TCGA). By projecting co-expression modules derived from bulk RNA-seq patient cohorts into single-cell datasets, we can layer disease-related information onto the single-cell dataset and attribute cell-state-specific expression patterns to the bulk RNA-seq data. We demonstrate projecting modules in this manner using co-expression modules from two bulk RNA-seq studies of AD, as the references and a human AD snRNA-seq dataset (57,950 nuclei from AD 11 samples and seven control samples of the PFC) as the query. These studies both used AD samples and controls from the same patient cohorts (Religious Orders Study and Memory and Aging Project, Mayo Clinic, Mount Sinai School of Medicine),,, but they took unique approaches for co-expression network analysis. The AMP-AD study from Wan et al. performed network analysis separately from each brain region, while, in our previous study, we performed consensus network analysis across the different brain regions. We projected these modules into a snRNA-seq dataset of AD and control samples from the PFC (Figure 7A), and we found distinct cell-type-specific expression patterns based on their MEs (Figures 7B and 7C). This analysis demonstrates hdWGCNA’s ability to transfer co-expression information across datasets to uncover otherwise unseen biological insights.

27Classical bioinformatic approaches for transcriptomics analysis such as differential gene expression are useful for finding individual genes that are altered in a particular disease or condition of interest, but they do not provide information about the broader context of these genes in specific pathways or regulatory regimes. For example, biological processes such as development or regeneration require coordination of distinct sets of genes in certain cell types with spatial specificity. Therefore, to understand these complex processes, we must look beyond individual genes. We developed hdWGCNA to provide a succinct methodology for investigating systems-level changes in the transcriptome in single-cell or ST datasets. We designed hdWGCNA to be highly modular, allowing for multi-scale analyses of different cellular or spatial hierarchies in a technology-agnostic manner.

28In this study, we demonstrated that hdWGCNA is compatible with single-cell and ST datasets and can be easily adapted for novel transcriptomics approaches such as ScISOrSeq. Co-expression networks have been successful for analyzing bulk proteomics datasets in human disease samples,, and we expect that hdWGCNA could be swiftly adapted for single-cell and spatial proteomics datasets as the technology matures and becomes more widely available. hdWGCNA includes built-in functions to leverage external biological knowledge sources to provide insight for co-expression networks, for example by comparing gene modules with functional gene sets such as disease-associated genes from GWAS expression quantitative trait loci (eQTLs), or transcription factor target genes. Unlike other network analysis pipelines such as single-cell regulatory network inference and clustering (SCENIC) or CellChat, hdWGCNA is a purely unsupervised approach and does not require prior knowledge or databases in the inference procedure. The co-expression information computed by hdWGCNA can be easily retrieved from the Seurat object to facilitate custom downstream analyses beyond the hdWGCNA package. hdWGCNA allows for comparisons between experimental groups via DME testing and module preservation analysis, which allowed us to identity inhibitory neuron modules that were dysregulated in ASD and enriched for ASD genetic risk genes, and microglial modules that were dysregulated in AD and enriched for DAM genes. Our network analyses of the ASD and AD datasets shows that hdWGCNA is capable of uncovering expanded disease-relevant gene sets via the interaction partners of known disease-associated genes such as the ASD SFARI genes or the AD DAM genes. We showed that the co-expression networks inferred by hdWGCNA were highly reproducible in unseen datasets, indicating that this is a robust methodology that reflects the underlying biology of the system of interest rather than picking up on technical artifacts. Further, hdWGCNA sheds new light on previously identified co-expression networks and gene modules by allowing modules to be projected from a reference dataset to a query dataset. The hdWGCNA R package directly extends the familiar Seurat pipeline and the SeuratObject data structure, enabling researchers to rapidly incorporate network analysis into their own workflows, going beyond cell clustering and differential gene expression analysis toward systems-level insights.

29Transcriptomic measurements of single cells are generally noisy, imposing challenges and limitations in the analysis of these datasets. Technical noise may arise from dropout events or from various steps in the experimental protocols, potentially making downstream data analysis and interpretation more difficult. hdWGCNA explicitly tries to handle the issues of technical dropouts and data sparsity by constructing networks in metacell or metaspot transcriptomic profiles rather than directly using the single-cell data. Furthermore, we show that module preservation statistical testing can assess the reproducibility of a co-expression network in external validation datasets, giving additional confidence in the results from hdWGCNA.

30Requests for further information should be directed to the lead contact, Vivek Swarup (vswarup@uci.edu).

31This study did not generate new unique reagents.

32Single-cell gene expression datasets typically contain many more zero valued entries than non-zero valued entries, meaning that these datasets are sparse. We formally define the sparsity of a gene expression matrix in Equation 1. Given an un-normalized counts matrix with genes and cells, sparsity is the sum of all zero valued elements.

33Complementing sparsity, the density of a single gene expression matrix is the sum of all non-zero valued elements divided by the total number of matrix elements, such that. A matrix is considered sparse if. Conventional single-cell gene expression assays yield sparse gene expression matrices. In general, correlations of sparse vectors may lead to downstream conclusions that are not robust or reproducible. Thus, as part of the hdWGCNA workflow, we propose a bootstrapped aggregation (bagging) algorithm to construct a gene expression matrix with considerably reduced sparsity prior to performing co-expression network analysis. Zero valued entries in a gene expression matrix have both biological and technical origins, and it is important to prioritize preserving relevant biological signals while reducing technical noise. For example, a biological zero may be attributed to a gene that is only expressed in a given cell population, whereas a technical zero may arise from low sequencing depth.

34We define the set of unique cell barcodes and the set of unique genes such that and. Transcriptomically similar cells are identified in a dimensionally-reduced representation of the gene expression matrix using the -nearest neighbors (KNN) algorithm, yielding sets of cells. Inherently, there is overlap between these sets of neighboring cells, and we include a parameter to control for the maximum allowable overlap. Cells are uniformly randomly sampled from, and gene expression signatures from are aggregated (sum or average) with their nearest neighbors. A cell is skipped if its neighbors have too much overlap with the set of neighbors from previously selected cells, in order to reduce redundancy in the downstream metacell expression matrix. The cell sampling loop converges when there are no more cells that satisfy the, or when the number of target metacells has been reached, yielding a metacell gene expression matrix. Sparsity of the input and output matrices and are computed to check that sparsity is reduced throughout this process. This metacell bagging algorithm is implemented as part of the hdWGCNA R package in the ConstructMetacells function, and the pseudocode for this algorithm is defined in Algorithm 1. We denote a vector containing the elements of the -th row of a matrix as and a vector containing the elements of the -th column as.Require: such that ⊳ gene expression matrix of genes and cellsRequire: such that ⊳ dimensional reduction of, with cells and dimensionsRequire: ⊳ the set of unique cell barcodesRequire:Require:Require: ⊳ is a matrix of rows and columns with the nearest neighbors of each cell ⊳ list containing barcodes of cells selected for aggregation, initialized as emptywhile and do ⊳ is randomly sampled from ⊳ the maximum number of overlapping neighbors between and barcodes inif.thenend if.end while. ⊳ subset of with the selected cells. ⊳ final metacell expression matrix

35Require: such that ⊳ gene expression matrix of genes and cells

36Require: such that ⊳ dimensional reduction of, with cells and dimensions

37Require: ⊳ the set of unique cell barcodes

38Require:

39Require:

40Require:

41⊳ is a matrix of rows and columns with the nearest neighbors of each cell

42⊳ list containing barcodes of cells selected for aggregation, initialized as empty

43while and do

44⊳ is randomly sampled from

45⊳ the maximum number of overlapping neighbors between

46and barcodes in

47if.

48then

49end if.

50end while.

51⊳ subset of with the selected cells

52. ⊳ final metacell expression matrix

53Sequencing-based ST approaches such as the Genomics Visium platform also yield sparse transcriptomic profiles, thus introducing the same potential pitfalls as single-cell data for co-expression network analysis. To alleviate these issues, we sought to develop a data aggregation approach similar to our metacell algorithm. This approach leverages spatial coordinates rather than the dimensionality-reduced representation (Figure S6). For each ST spot, we obtain a list of physically neighboring spots. We then devise a grid of ”principal spots”, which are evenly spaced spots throughout the input tissue which serve as anchor points for aggregating neighboring spots. Each principal spot and its neighbors are aggregated into one metaspot, with at most seven spots merging into one metaspot and at most two overlapping spots between metaspots. We implemented this procedure as part of the hdWGCNA R package in the MetaspotsByGroups function. Similar to the MetacellsByGroups function, the user may specify groups within the Seurat object to perform the aggregation, such that metacells would only be grouped within the same tissue slice, anatomical region, or other annotation. For all downstream analysis with hdWGCNA, the metaspot expression dataset can be used in place of the metacell expression matrix.

54Following metacell or metaspot construction, hdWGCNA constructs co-expression networks and identifies gene modules, building off of the WGCNA workflow.,,, The gene-gene adjacency matrix is computed by taking the pairwise correlation of genes in in the metacell expression matrix, or in a subset of for a specified cell population. Consider the gene expression vectors and for an arbitrary pair of genes, we compute the signed correlation as:

55Note that is a linear transformation that retains the sign of the correlation while satisfying. We define as a symmetric adjacency matrix of size containing the signed correlations for all pairs as in Equation 2. In order to emphasize strong correlations, we raise the elements of to a power, and we refer to this as soft power thresholding.

56Now we have the gene-gene correlation raised to a power, and an alternative metric which also retains the sign of the correlation between these genes. The final co-expression network is then computed as a signed topological overlap matrix (TOM). The TOM describes shared neighbors between the a pair of genes. We define the signed TOM aswhere and represent the connectivity between genes and In the signed TOM, negative correlations serve to negatively reinforce the network connection, which is not the case in the unsigned TOM.

57Genes are then grouped into modules based on the TOM network representation using the Dynamic Tree Cut algorithm, such that co-expression modules consist of genes with high topological overlap. Dynamic Tree Cut hierarchically clusters genes based on their dissimilarity in the TOM, denoted as, thereby yielding a mapping between module assignments and gene names. The overall process transforming a metacell expression matrix to a signed TOM co-expression network is implemented as part of the hdWGCNA R package in the ConstructNetwork function. Here we described the recommended workflow, using a signed adjacency matrix and a signed TOM, but ConstructNetwork can optionally construct unsigned or signed hybrid networks as well.

58Module eigengenes (MEs) are a convenient metric to summarize the gene expression of a given co-expression module. While the co-expression network was computed using the metacell expression matrix, we compute MEs in the single-cell expression matrix, thus yielding information about the activity of each module in each cell. The expression matrix for the -th module consisting of genes is. The ME for module is then computed by performing singular value decomposition (SVD), such that. Prior to running SVD, must be scaled and centered, and we accomplish this using the Seurat function ScaleData. Importantly, ScaleData enables us to optionally perform regression to diminish the effects of selected technical covariates prior to computing MEs. The first column of, containing the right-singular vectors, is the ME of module.

59While SVD or other dimensionality reductions on a single-cell gene expression matrix contains critical biological information, technical artifacts are also present in these representations. There are many computational methods aiming to reduce technical effects in a reduced dimensional space, and these methods are often referred to as “batch-correction” or “integration” approaches. In particular, Harmony is an algorithm well suited for correcting batch effects that may be present in a dimensionality-reduced single-cell expression dataset, and here we propose applying Harmony to MEs to maximize the biological information content of each ME. We implemented the ME computation algorithm, as defined in Algorithm 2, as part of the hdWGCNA R package in the function ModuleEigengenes.Require: such that ⊳ normalized gene expression matrix of genes and cellsRequire: modules ⊳ the table containing mappings between genes and modules.Require: mods ⊳ list of modules.Require: covariates ⊳ covariates to regress.Require: batches ⊳ batch identity to correct with Harmony, or null to ignore.ME fordomodules subset(modules, module = = ) modules [,gene] ScaleData(, covariates) SVD ME ifthen Harmony(, batches)ME end if.ME [ ME, ME ]end for.

60Require: such that ⊳ normalized gene expression matrix of genes and cells

61Require: modules ⊳ the table containing mappings between genes and modules.

62Require: mods ⊳ list of modules.

63Require: covariates ⊳ covariates to regress.

64Require: batches ⊳ batch identity to correct with Harmony, or null to ignore.

65ME

66fordo

67modules subset(modules, module = = )

68modules [,gene]

69ScaleData(, covariates)

70SVD

71ME

72ifthen

73Harmony(, batches)

74ME

75end if.

76ME [ ME, ME ]

77end for.

78In a typical hdWGCNA workflow, we perform metacell bagging, co-expression network analysis, module identification, and ME computation using the same single-cell gene expression dataset, starting from the expression matrix. Given the module-gene assignment table derived from a reference dataset, we can run the ModuleEigengenes algorithm on a query dataset where the genes in must be contained in the set of genes in such that. We implemented this process in the hdWGCNA R package as the ProjectModules function. Importantly, we designed ProjectModules to be agnostic towards the data modality or species used in the reference and query datasets, thereby allowing for a host of comparative analyses. ProjectModules can facilitate cross-species analysis leveraging a table that maps gene symbols between two genomes. Modules can be projected into epigenomic data modalities such as single-cell assay for transposase accessible chromatin with sequencing (scATAC-seq) provided a measure of gene expression estimated from chromatin accessibility, such as Signac gene activity or ArchR gene scores. This approach can also be used to project modules from bulk expression datasets into single-cell or spatial transcriptomics datasets.

79hdWGCNA greatly extends upon scWGCNA, our previous method for co-expression network analysis in single-cell transcriptomics data. scWGCNA was originally used to identify co-expression networks using bulk and single-cell RNA-seq together, and in another study we showed that scWGCNA was suitable for network analysis using scRNA-seq alone. Contrasting the hdWGCNA package, the implementation of scWGCNA was an R package containing a single function for metacell construction, and a single tutorial to cover the basics of network analysis using the WGCNA package with the metacell matrix. We implemented hdWGCNA as an open-source object-oriented R package that leverages the widely used SeuratObject data structure. The hdWGCNA R package includes all necessary functions for network inference, data visualization, statistical testing, and downstream analysis such as pathway enrichment. Further, hdWGCNA includes functions to extract the network data from the SeuratObject to easily facilitate custom analysis with external Bioconductor or R packages. In order for hdWGCNA to be widely useful across the genomics community, we developed a detailed documentation website containing tutorials for network analysis in single-cell and spatial transcriptomics data, as well as tutorials for advanced analysis like consensus network analysis and network preservation testing. Unlike scWGCNA, the metacell construction algorithm (Algorithm 1) in hdWGCNA includes new parameters to avoid redundant metacells, the module eigengene algorithm in hdWGCNA (Algorithm 2) accounts for batch effects and additional covariates in the input dataset, and hdWGCNA contains functions to handle spatial transcriptomics datasets. Several key steps in co-expression network analysis, like calculating module eigengenes and eigengene-based connectivity, have been re-implemented to operate on sparse matrices, greatly decreasing runtime and memory usage. hdWGNCA is completely technology agnostic, and can be adapted to handle high dimensional counts matrices from any single-cell or spatial transcriptomics platform. Additionally, hdWGCNA includes a novel approach for visualizing genes and the underlying network in a two-dimensional manifold of co-expression space using UMAP. As shown throughout this manuscript, hdWGCNA includes functions for projecting co-expression networks into a variety of external datasets. The widespread adoption of single-cell genomics has led to many biologists running their own computational analysis, and we designed the hdWGCNA R package with these individuals in mind through our various step-by-step tutorials and detailed documentation.

80The key resources table details the different datasets used throughout this manuscript. We used several published datasets generated by our own group,,, and sequencing data was not re-downloaded for these studies. For all human snRNA-seq datasets, we applied a uniform processing pipeline to process each dataset starting from the raw sequencing data and resulting in an anndata object containing UMI counts, normalized gene expression, cluster identities, and cell type annotations. Parameters used throughout this processing pipeline vary slightly between different datasets, and all parameters are noted in the data processing scripts in our github repository. For each biological replicate, we used the kb count function from kallisto bustools to psuedoalign raw sequencing reads to the reference transcriptome and quantify gene expression attributed to each cell barcode. The human reference transcriptome (GRch38) was obtained from the Genomics website (version 2020-A, July 2020), and was re-formatted for use with kallisto bustools using the kb ref function. For each of the UMI counts matrices, we used the remove-background function from cellbender to simultaneously identify which barcodes corresponded to cells and to remove counts attributed to ambient RNA. We then used scrublet to compute ”doublet scores”, the likelihood of each barcode mapping to more than one cell. Counts matrices from each biological replicate in a given dataset are then merged into a single anndata object, and any relevant sample level meta-data (age, sex, disease status) was stored in the adata.obs table. We performed a percentile filtering of cells that were outliers from each dataset based on the number of UMI per cell the percentage of UMI attributed to mitochondrial genes per cell, and the doublet score. Filtering based on these criteria was performed in each sample, as well as dataset-wide. After filtering, downstream data processing steps were carried out with SCANPY. The UMI counts matrix was normalized with using the functions sc.pp.normalize_total and sc.pp.log1p. Highly variable genes were identified using the function sc.pp.highly_variable_genes, and these genes are used as the features for downstream analysis steps such as principal component analysis (PCA). The normalized expression matrix was then scaled to unit variance and centered at zero using the function sc.pp.scale. PCA was performed on the scaled expression matrix using the function sc.tl.pca. Harmony was used to correct the PCA matrix for batch effects using the function sc.external.pp.harmony_integrate. The harmonized PCA matrix was then used to construct a cell neighborhood graph using the function sc.pp.neighbors. The cell neighborhood graph was then used to compute a two-dimensional representation of the data with uniform manifold approximation and projection using the function sc.tl.umap, and to group cells into clusters with Leiden clustering using the function sc.tl.leiden. We inspected the gene expression signatures in each Leiden cluster for a panel of canonical cell-type marker genes in order to assign a cell-type label to each cluster, and to identify additional doublet clusters that may have escaped the previous filtering steps. The distribution of quality control metrics was inspected in each cluster. We filtered out cells belonging to clusters that displayed conflicting expression of cell-type marker genes, or were outliers in their quality control metrics. After filtering these low-quality clusters, we ran UMAP and Leiden clustering again, resulting in the final processed dataset. We used a custom script to convert the datasets from anndata to SeuratObject by saving the individual components (counts matrix, cell meta-data, gene meta-data, dimensionality reductions, etc.) in Python and then loading them back into R to create a SeuratObject.

81We performed an iterative co-expression network analysis of the major cell types (ASC, EX, INH, MG, ODC, OPC) in the human PFC snRNA-seq dataset from Zhou et al., only including samples from control brains (36,671 cells and 36,601 genes). We retained genes that were expressed in at least 5% of cells for downstream analysis. Metacells were computed separately for each major cell type and each sample using the hdWGCNA function MetacellsByGroups, aggregating 25 cells per metacell. Further, we ran MetacellsByGroups while varying the parameter in order to asses the resulting metacell expression matrix sparsity. For each cell type, we applied the following hdWGCNA commands with default arguments to perform network analysis: TestSoftPowers, ConstructNetwork, ModuleEigengenes, ModuleConnectivity, and RunModuleUMAP. We performed module preservation analysis of the ODC co-expression modules in an external snRNA-seq dataset of the human PFC. Modules were projected from the reference to query dataset using the hdWGCNA function ProjectModules, and the module preservation test was performed using ModulePreservation with 100 permutations.

82For the purpose of co-expression network analysis, we compared our metacell aggregation approach (Algorithm 1) with two alternative approaches, namely Metacell2 and SEACells. We ran the three metacell approaches using the recommended settings on the same dataset, and then ran hdWGCNA on each of the resulting metacell expression matrices. We used a scRNA-seq of 6,800 CD34+ hematopoietic stem and progenitor stem cells included with the SEACells package, and we used the cluster annotations from the original study. Notably, SEACells and Metacell2 do not account for cell labels in their aggregation procedures, which may result in a number of metacells containing transcriptomes from differently labeled cells. For the hdWGCNA metacell algorithm, we aggregated 50 cells per metacell. For the three metacell expression matrices derived from the different algorithms, we performed co-expression network analysis with the standard hdWGCNA pipeline by sequentially running the following functions with default parameters: TestSoftPowers, ConstructNetwork, ModuleEigengenes, ModuleConnectivity, and RunModuleUMAP. With the same cluster settings, Dynamic Tree Cut recovered a different number of co-expression modules for the three methods (hdWGCNA: 16 modules; MC2: 13 modules; SEACells: 20 modules). We performed pairwise comparisons between the gene modules detected with each metacell approach using Fisher’s exact test to test module overlaps. Additionally, we performed rank-rank hypergeometric overlap (RRHO) tests using the RRHO function from the R package RRHO (version 1.13.0) to compare the kME ranking between modules across methods. To compare MEs and Seurat module scores, we ran the AddModuleScore function, and computed Pearson correlations between each ME and each module score.

83We obtained a publicly available scRNA-seq dataset from Parse Biosciences of 1M peripheral blood mononuclear cells (PBMCs) from twelve healthy donors and twelve Type-1 diabetic donors generated using the Evercode Whole Transcriptome Mega protocol. This analysis was performed on a compute cluster with 200 GB of memory and eight CPU cores. The UMI counts matrix and sample meta data was downloaded from Parse Biosciences’ Website. We processed the counts matrix using SCANPY using a similar pipeline as described in the reprocessing published dataset section. For quality control, we excluded cells with greater than 25% mitochondrial reads, greater than 5,000 genes, and greater than 25,000 counts. After dimensionality reduction with PCA, Harmony batch correction, and Leiden clustering (resolution = 1), we annotated cell populations using PBMC marker genes obtained from Azimuth. We excluded clusters with conflicting cell-type markers as potential doublet populations, retaining a total of 965,363 cells and 26,862 genes for downstream analysis. The major cell compartments recovered in this analysis were similar to those reported by Parse Biosciences in their analysis, including as T-cells, B-cells, monocytes, dendritic cells, basophils, and plasmablasts. Following the SCANPY data processing, we wrote the individual components (counts matrix, cell meta-data, gene meta-data, dimensionality reductions, etc.) to disk so they could be loaded into R and assembled into a Seurat object.

84We performed co-expression network analysis iteratively for the plasmablast, T-cell, B-cell, monocyte, and dendritic cell compartments using an hdWGCNA pipeline for each group (Figure S5). Metacells were constructed separately for each sample and each cell cluster with the hdWGCNA function MetacellsByGroups, aggregating 50 cells per metacell. The metacell aggregation step had a runtime of 85 min and 59 s. For each cell population, we first subset the Seurat object for the cell population of interest and then performed the standard hdWGCNA pipeline by sequentially running the following functions with default parameters: TestSoftPowers, ConstructNetwork, ModuleEigengenes, ModuleConnectivity, and RunModuleUMAP. We note that for the largest cell population (T-cells, 555,417 cells), the runtime for the network construction step was 186 s.

85We tested the runtime and memory usage of the primary co-expression network analysis functions in hdWGCNA using the Velmeshev et al. 2019 dataset. We selected the neuronal cell population from the dataset for network analysis, and downsampled the dataset at different sizes ranging from 1,000 to 50,000 cells to test the runtime and memory usage as a function of the number of cells in the input dataset. The following functions were tested: SetupForWGCNA, MetacellsByGroups, TestSoftPowers, ModuleEigengenes, and ModuleConnectivity. We tested ModuleEigengenes with and without Harmony correction. All of these tests were done using eight parallel threads, and hdWGCNA can be sped up further by increasing the number of parallel threads. Importantly, the number of input genes and other network analysis parameters also have an effect on runtime and memory usage.

86We tested the functional coherence of hdWGCNA co-expression networks using the Extending ’Guilt-by-Association’ by Degree (EGAD) algorithm. Connected genes in biological networks are potentially involved in the same processes, and EGAD evaluates this network property given a set of gene-process annotations. We performed functional coherence testing with the EGAD R package (version 1.18.0) using the six cell-type-specific co-expression networks from the Zhou et al. 2020 human PFC dataset. We downloaded a table of gene ontology associations for each gene from ensembl biomart, and formatted this table using the EGAD function make_annotations. We then ran the functional coherence test with EGAD using the function run_GBA, using the TOM as the input network, and we report the distributions of area under the receiver operating characteristic curve (AUC) values for each tested biological process in the six co-expression networks.

87We used the xgboost R package (version 1.7.3.1) to perform XGBoost regularized regression analysis to predict a given gene’s expression based on the expression of the top ten module hub genes for the module each gene was assigned to. This analysis was done using the six cell-type-specific co-expression networks described in the iterative network analysis of major cell types in the human cortex section. We performed 5-fold cross validation, and measured the performance of the model as a test set root-mean-square error (RMSE) averaged across the 5-folds. We ran XGBoost for 100 iterations for each individual test, with a maximum tree depth of 3 and regularization alpha of 0.5.

88We collected the publicly available Genomics Visium mouse brain dataset using the SeuratData R package. This dataset consists of an anterior and a posterior slice from a sagittal brain section, which we merged into a single Seurat object comprising 6,049 ST spots and 31,053 genes. We processed this dataset using the standard Seurat pipeline by sequentially running the following commands: NormalizeData, FindVariableFeatures, ScaleData, RunPCA, FindNeighbors, FindClusters, and RunUMAP. The top thirty PCs were used for Louvain clustering and UMAP. While ST spots were clustered based on transcriptomic information alone, we were able to annotate them based on anatomical features.

89Neighboring ST spots were aggregated into metaspots in the anterior and posterior slices using the hdWGCNA function MetaspotsByGroups. We retained genes expressed in 5% of spots for downstream analysis, totaling 12,355 genes. We tested for the optimal soft-power threshold based on the the fit to a scale-free topology using the hdWGCNA function TestSoftPowers. The co-expression network was constructed using all ST spots spanning both the anterior and posterior slices using the hdWGCNA function ConstructNetwork with the following parameters: networkType = “signed”, TOMType = “signed”, soft_power = 5, deepSplit = 4, detectCutHeight = 0.995, minModuleSize = 50, mergeCutHeight = 0.2. Module eigengenes and eigengene-based connectivities were computed using the ModuleEigengenes and ModuleConnectivity functions respectively. This approach identified 12 spatial co-expression modules, and we visualized the spatial distributions of these modules by plotting their MEs directly onto the biological coordinates for each spot. The co-expression network was projected into two dimensions using UMAP with the hdWGCNA function RunModuleUMAP, and we used the top five hub genes (ranked by kMEs) as the input features for UMAP. We used the R package enrichR (version 3.0) to perform enrichment analysis on the top 100 genes in each module ranked by kME using the following databases: GO_Biological_Process_2021, GO_Cellular_Component_2021, GO_Molecular_Function_2021, WikiPathway_2021_Mouse, and KEGG_2021_Mouse. We assessed the overlap between genes from these spatial co-expression modules and differentially expressed genes in each cluster from a recent snRNA-seq study of the whole mouse brain using Fisher’s exact test implemented in the R package GeneOverlap (version 1.26.0). Finally, we performed a separate network analysis on a subset of the ST dataset only containing the cortical layers 2–6, and we followed an identical hdWGCNA analysis pipeline to the full ST dataset for the cortical analysis.

90We performed isoform co-expression network analysis in radial glia lineage cells (radial glia, astrocytes, ependymal cells, and neural intermediate progenitor cells) from mouse hippocampus ScISOrSeq dataset from Joglekar et al. using the hdWGCNA R package. The gene-level counts matrix for this dataset was obtained from the Gene Expression Omnibus database (GEO: GSE15845), and the isoform-level counts matrix was obtained directly from the authors of the original study. We formatted this dataset as a Seurat object with an isoform-level expression assay and a gene-level expression assay. The standard Seurat processing pipeline was used on the gene-level expression assay, where we sequentially ran the functions NormalizeData, FindVariableFeatures, ScaleData, and RunPCA with default parameters. The dataset was projected into two dimensions by running UMAP on the PCA matrix with 30 components using the RunUMAP function. For all downstream purposes, the cell-type annotations from the original study were used.

91Radial glia cells were selected for network analysis, and isoforms expressed in fewer than 1% of these cells were excluded, yielding a set of 2,190 cells and 10,375 isoforms from 4,770 genes. We constructed metacells separately for each cell type on the isoform-level expression assay using the hdWGCNA function MetacellsByGroups with. We performed a parameter sweep for the soft-power threshold using the function TestSoftPowers. The isoform co-expression network was constructed using the ConstructNetwork function with the following parameters: networkType = “signed”, TOMType = “signed”, soft_power = 5, deepSplit = 4, detectCutHeight = 0.995, minModuleSize = 50, mergeCutHeight = 0.5. This approach identified 11 isoform co-expression modules. Isoform-level module eigenisoforms were computed using the ModuleEigengenes function, and eigenisoform-based connectivity was computed using the ModuleConnectivity function with default parameters. We computed a semi-supervised UMAP projection of the co-expression network using the hdWGCNA function RunModuleUMAP, with the module labels and the top six hub isoforms (by kMEiso) per module as the input features. We used the enrichR to identify enriched pathways in each module ranked by using the following databases: GO_Biological_Process_2021, GO_Cellular_Component_2021, GO_Molecular_Function_2021, WikiPathway_2021_Mouse, and KEGG_2021_Mouse.

92To assess isoform co-expression network dynamics throughout the cellular trajectories within the radial glia lineage, we performed pseudotime analysis using Monocle 3 (version 1.0.0). We computed a UMAP of just radial glia lineage cells using the Monocle 3 function run_umap. A trajectory graph was built on this UMAP representation using the function learn_graph, and pseudotime was calculated with the function order_cells using the radial glia cells as the starting point. We split the pseudotime trajectory into three lineages based on the distinct cell fates (astrocyte, neuronal, and ependymal). We grouped cells into 50 evenly-sized bins throughout each trajectory, and we applied loess regression to the average module eigenisoform of each module in these bins to inspect the dynamics of each module throughout development. We wrote a custom script to generate a GTF of isoform models output from the ScISOrSeq pipeline. To visualize expressed isoforms, we plotted isoforms from this GTF on the UCSC genome browser as well as in Swan.

93We selected inhibitory neurons from the Velmeshev et al. human autism spectrum disorder (ASD) snRNA-seq dataset for co-expression network analysis. Of the 121,451 cells in this dataset, 20,249 were labeled as inhibitory neurons based on marker gene expression profiles. We retained 11,194 genes which were expressed in at least 10% of cells from any cluster, and had non-zero variance in the inhibitory neuron population. Metacell transcriptomic profiles were constructed separately for each of the 54 samples and each cell type using the hdWGCNA function MetacellsByGroups, aggregating 50 cells into one metacell. We selected a soft-power threshold based on the parameter sweep performed with the TestSoftPowers function. The co-expression network was computed with the ConstructNetwork function with the following parameters: networkType = “signed”, TOMType = “signed”, soft_power = 9, deepSplit = 4, detectCutHeight = 0.995, minModuleSize = 50, mergeCutHeight = 0.2. Module eigengenes were computed using the ModuleEigengenes function, and we applied Harmony to correct MEs based on sequencing batch. Eigengene-based connectivity for each gene was computed using ModuleConnectivity. The co-expression network was embedded in two dimensions using UMAP with the RunModuleUMAP function with the top five genes (ranked by kMEs) per module as the input features. Distributions of MEs were compared between ASD and control samples for each inhibitory neuron subpopulation using a two-sided Wilcoxon rank-sum test with the R function wilcox.test. We used the enrichR to perform enrichment analysis on the top 100 genes in each module ranked by kME using the following databases: GO_Biological_Process_2021, GO_Cellular_Component_2021, GO_Molecular_Function_2021, WikiPathway_2021_Human, and KEGG_2021_Human. Furthermore, we computed the overlap between co-expression modules and ASD-associated genes from the SFARI Gene database using the R package GeneOverlap, which calculates the overlap between sets of genes using Fisher’s exact test.

94We performed consensus co-expression network analysis of microglia in Alzheimer’s disease (AD) using three published snRNA-seq datasets.,, The individually processed datasets were merged into a single Seurat object comprising 189,127 nuclei, and the datasets were integrated into a common dimensionally-reduced space using PCA and Harmony. We retained all nuclei labeled microglia for network analysis based on expression of canonical marker genes such as CSF1R (9,904 nuclei), and genes expressed in at least 5% of microglia from any of the three studies were retained (7,900 genes). Metacells were constructed in groups of cells based on AD diagnosis status and study of origin, aggregating 25 cells per metacell. Within hdWGCNA, we used the SetMultiExpr function to create a list of expression matrices containing the selected genes and metacells for the three studies. We performed a separate parameter sweep for the three expression matrices using the hdWGCNA function TestSoftPowerConsensus, ensuring that we used an appropriate value for each dataset (Mathys et al.:, Zhou et al.:, Morabito & Miyoshi et al.: ). The consensus co-expression network was contructed using the hdWGCNA function ConstructNetwork using the consensus = TRUE option. Individual TOMs were computed for each dataset, and they were scaled based on the 80th percentile in order to alleviate different statistical properties specific to each dataset rather than the underlying biology. A consensus TOM was computed by taking the element-wise minimum of the individual TOMs from each dataset. Therefore, large topological overlap values between two genes, which indicate a strong co-expression relationship, are supported across all three datasets in the consensus TOM. We performed hierarchical clustering on the consensus TOM, and we used the Dynamic Tree Cut algorithm was used to identify consensus co-expression modules based on the hierarchy. Module eigengenes were computed using the ModuleEigengenes function, and we applied Harmony to correct MEs based on the dataset of origin. Eigengene-based connectivity for each gene was computed using ModuleConnectivity. We visualized the network using UMAP with the top ten hub genes (ranked by kMEs) per module as the input features, annotating the hub genes and known disease-associated microglia genes. We used the enrichR to perform enrichment analysis on the top 100 genes in each module ranked by kME using the following databases: GO_Biological_Process_2021, GO_Cellular_Component_2021, GO_Molecular_Function_2021, WikiPathway_2021_Human, and KEGG_2021_Human.

95We sought to model the transcriptional dynamics governing the shift between homeostatic and activated microglia in AD, therefore we performed pseudotime analysis using Monocle 3 to build a continuous trajectory of microglia cell states. A trajectory graph was built on the microglia UMAP using the function learn_graph, and pseudotime was calculated with the function order_cells. We oriented the start of pseudotime based on the expression of homeostatic microglia marker genes, such as P2RY12, CX3CR1, and CSF1R. We grouped cells into 50 evenly-sized bins throughout each trajectory, and we applied loess regression to the average module eigengene of each module in these bins to inspect the dynamics of each module throughout the microglia trajectory.

96To link the integrated microglia snRNA-seq dataset with polygenic risk of disease for individual cells, we used the scDRS python package (version 1.0.0). This pipeline takes 1) a set of putative disease genes derived from GWAS summary statistics and 2) a scRNA-seq dataset as inputs, and outputs disease enrichment statistics for a given disease (raw and normalized disease scores, cell-level scDRS p value, and -scores converted from the p values). GWAS summary statistics of 74 diseases and complex traits supplied by scDRS were utilized as gene sets, among which a gene set by Jansen et al. provided the set of genes associated with AD. We then visualized the AD scDRS -scores in the integrated AD microglia trajectory, and we correlated the scDRS score with the trajectory using a Pearson correlation.

97We performed module preservation analysis in a variety of external datasets from human and mouse,,,,, to test for the reproducibility of the consensus AD microglia modules in the microglia population from each dataset. We used the hdWGCNA function ProjectModules to compute module eigengenes for the consensus AD microglia modules for each query dataset. The module preservation test was performed using the hdWGCNA function ModulePreservation with 100 permutations, and we reported the preservation -summary statistics in a heatmap. For the Morabito & Miyoshi et al. snATAC-seq dataset, we used the gene activity representation as a gene-level summary of chromatin accessibility in order to assess the module preservation at the epigenomic level.

98We projected gene co-expression modules from two bulk RNA-seq studies of AD, into a published snRNA-seq study of AD to assess their expression patterns within various cell populations. While both of these studies used the samples from the same bulk RNA-seq cohort, the set of modules from Morabito et al. 2020 was based on a consensus network analysis across six brain regions while the other set of modules from the AMP-AD study were constructed separately for seven different brain regions. Module eigengenes were computed for each of these bulk RNA-seq modules in the snRNA-seq dataset using the hdWGCNA function ProjectModules, using Harmony to correct MEs based on sequencing batch. We visualized the MEs of the projected modules in the snRNA-seq dataset using the Seurat function DotPlot.