HCNetlas: A reference database of human cell type-specific gene networks to aid disease genetic analyses

We employed both single-cell and single-nucleus RNA sequencing (scRNA-seq and snRNA-seq) data from multiple sources to reconstruct HCNetlas CGNs. We acquired scRNA-seq data for 329,762 immune cells spanning 16 tissues from 12 deceased donors in a cross-tissue immune cell atlas [10]. We used the author-provided cell type labels, which were initially annotated by CellTypist and subsequently underwent manual curation. Additionally, we integrated data from Tabula Sapiens [11], which included the human transcriptome reference for 249,961 immune cells across 24 tissues from 15 donors. For constructing brain CGNs, we utilized snRNA-seq data of Allen Brain Atlas (http://www.brain-map.org) [12], derived from 76,533 nuclei in the primary motor cortex (M1) and 166,868 nuclei in the middle temporal gyrus (MTG) using 10× Genomics Chromium platform. All data were processed through alignment, quantification by Cell Ranger, and cell type annotation.

We harmonized pre-annotated cell types across all data sets and applied scHumanNet (v. 1.0.0) [5] to the single-cell transcriptomic data to generate CGNs for various tissues and cell types. Briefly, scHumanNet takes a single-cell transcriptome count matrix for a specific cell type, followed by data imputation, transformation, and normalization. It then produces a transformed gene activity score matrix and filters the reference interactome, HumanNet, for gene pairs in the matrix that pass a minimum activity score, creating a CGN. The edge weight scores are inherited from the reference interactome. scHumanNet has been shown to be a superior network inference model compared to other single-cell based methods. To construct CGNs for major cell types, we aggregated sub-cell types into broader categories, such as T cells, B cells, and myeloid cells. In total, 198 CGNs were generated, encompassing 25 tissues and 61 cell types.

We performed dimension reduction analysis with uniform manifold approximation and projection (UMAP) to visualize interrelationships among HCNetlas CGNs. Utilizing scHumanNet for CGN construction, which references the HumanNet [6] comprising 18,593 human genes, we generated binary profile vectors indicating the presence (1) or absence (0) of each gene in the network for every CGN. For visualizations, we employed the UMAP package in R (v. 0.2.10), setting min_dist to 0.5 to balance the trade-off between local and global structure in the data.

To determine the cell type-specific functionality of the HCNetlas CGNs, we explored the enrichment of cell type-associated genes, particularly for B cells and T cells. We collated cell type-associated genes from 3 authoritative databases: Gene Ontology biological process (GOBP), CellMarker, and the Azimuth cell type database [13–15]. We postulated that genes functionally connected within a CGN reflect the properties of their respective cell type, thus we considered the interconnectivity within genes for each cell type as a measure of the cell type specificity.

We assessed whether HCNetlas CGNs can discern connections between diseases and cell types by profiling CGNs with disease-association scores. Within each CGN, we ranked the 18,593 genes from the HumanNet reference interactome by degree centrality, utilizing the GetCentrality() function of the scHumanNet. Disease gene sets, totaling 5,763, were compiled from DisGeNET [16] and GWAS Catalog [17] for the analysis. Assessment of disease gene set association with each CGN was performed with ssGSEA (v. 2.0) [18] and GSVA [19], the latter via the gsva() function from the GSVA package (v. 1.38.2). Hierarchical clustering of disease gene sets and CGNs was performed using the R package ComplexHeatmap, with the default parameters, where cluster.rows and cluster.cols were both set to TRUE.

We acquired scRNA-seq data for peripheral blood mononuclear cells (PBMCs) from 41 patients with SLE and 15 healthy controls as reported by Nehar-Belaid and colleagues [20]. To ensure a consistent analysis, we excluded 2 SLE patient samples lacking SLEDAI scores. Following quality control measures, including the removal of doublets using the DoubletFinder (v. 2.0.0) [21] package and the exclusion of cells with fewer than 400 transcripts or over 5% mitochondrial gene content, the data set was narrowed to approximately 276,000 cells. After normalization and scaling with Seurat package (v. 4.1.1), we identified 3,000 variable genes using the vst() and FindVariableFeatures() functions. Batch effects were mitigated by applying principal component analysis (PCA) and Harmony (v. 1.0) [22] (dims = 40), and cellular clustering was performed using the Louvain method (resolution = 1.5), followed by UMAP visualization with 40 dimensions. Cell types were manually annotated using canonical markers after optimizing the number of principal components and clustering resolution.

For constructing SLE-specific CGNs, we focused on the curated data from SLE patients. We built networks for B cells, CD4⁺ T cells, CD8⁺ T cells, myeloid cells, and NK cells. Using the Compactness() function, we performed differential compactness analysis. We referenced 184 SLE-associated genes from KEGG pathway (I05322) and KEGG disease (H00080) databases, comparing connectivity within disease CGNs and HCNetlas CGNs, and visualized the networks in Cytoscape (v 3.9.1) [23].

Node centrality within these networks was computed using the GetCentrality() function from the scHumanNet package. We compared the percentile ranks of centrality for disease CGNs against the reference CGNs of HCNetlas using the DiffPR.HCNetlas() function. Genes showing differential hubness were pinpointed with FindDiffHub.HCNetlas(), with significance defined by a q-value

We compiled interferon-stimulated genes (ISGs) from hallmark gene sets of the molecular signature database (MSigDB) and the Immunological Genome Project (ImmGen) [24], resulting in a total of 423 ISGs. The efficacy of prediction for ISGs by hubness within SLE-associated CGNs was assessed using receiver operating characteristic (ROC) curve analysis. The ROC curve was generated using the roc() function from the pROC package (v. 1.18.0).

In our study of AD, we used snRNA-seq data from 12 patients with annotated cell types from Morabito and colleagues [25]. Since the data were derived from the prefrontal cortex of brain tissues, the generated CGNs for AD were compared with reference CGNs for the primary motor cortex (M1) from the HCNetlas. We grouped the cell type annotations into 4 main categories: astrocytes, inhibitory neurons, excitatory neurons, and oligodendrocytes. The identification of differential hubness genes and DEGs within these cell types was carried out using the same methodology applied in the analyses of SLE. To ascertain the relevance of AD-associated genes predicted by our differential hubness analysis, we referenced genes linked to AD in the KEGG pathway (M16024), MSigDB (M35868), and Wightman and colleagues [26].

Considering the association of gain-of-hubness and loss-of-hubness genes with AD in inhibitory and excitatory neurons, we constructed network-ranked signatures for both reference and AD-specific CGNs for the cell types. The signature genes were based on the top 10 hub genes by degree centrality within each CGN. The networks of these top-tier hub genes were visualized using the Cytoscape software [23]. Furthermore, we conducted gene set enrichment analysis (GSEA) on these network-ranked signatures using the enrichR package [27]. To evaluate the pathways differentially associated between disease CGNs and reference CGNs, we introduced a metric called diffQ, calculated as follows:

To create lung cancer-specific CGNs, we used scRNA-seq data from 29 tumor tissues provided by Qian and colleagues [28]. We retained the pre-annotated cell type identifications from the data sets. For comparison with reference CGNs derived from paired normal tissues, we constructed networks from both the lung cancer and healthy control data. Following the scHumanNet protocol, we generated networks and defined differential hubness genes using FindDiffHub(). The process of identifying differential hubness genes within each cell type was conducted using the same methodology employed in the SLE analyses. Similarly, the identification of DEGs followed the methodology used in the SLE studies, with the exception that genes exhibiting an absolute log₂-fold change

To validate the lung cancer relevance of the identified genes, we referenced the Cancer Gene Census, CancerMine, and IntOGen databases [29–31]. We assessed the proportion of lung cancer-associated genes detected uniquely through differential hubness, uniquely through DEGs, and by the intersection of both methods. Furthermore, we analyzed 42 immune checkpoint molecules listed by Auslander and collagues [32] to determine if cell type-specific genes vital for cancer immunity are discernible through both expression-based and network-based analyses.

We investigated the prognostic potential of genes identified by cell type-specific differential hubness and differential expression using survival analysis on TCGA lung cancer data sets (TCGA-LUSC, TCGA-LUAD). Initially, we identified a total of 379 gain-of-hubness genes and 211 up-regulated DEGs from 3 major cell types: B cells, T cells, and myeloid cells. Subsequently, genomic and clinical data for 1,017 lung cancer samples were acquired from the GDC portal [33]. The STAR-Counts data underwent preprocessing, log-normalization, and variance stabilization using the vst() function in the DESeq2 R package (v. 1.30.1). With the application of GSVA [19], we evaluated the association of both gain-of-hubness genes and up-regulated DEGs with each tumor expression profile of the patients. Patients were then classified into upper and lower quartile groups based on their GSVA scores. These groups were further examined through Kaplan–Meier survival curves. To ensure the reliability of our findings, we adjusted all p-values obtained from the survival analysis using the Benjamini–Hochberg method to control the FDR.

To build reference CGNs, we utilized scRNA-seq data from immune cell atlas project [10] and The Tabula Sapiens [11], as well as single-nucleus RNA-sequencing (snRNA-seq) data from the Allen Brain Atlas [12]. Our single-cell transcriptomic data set comprised 763,559 cells from 28 donors. We generated gene networks for each predefined cell type using the scHumanNet framework [5] (Fig 1A), providing a comprehensive baseline for identifying disease-associated genes and cell types.

Fig 1. Overview of HCNetlas.

(A) Schematic representation of the workflow from single-cell transcriptomic data collection to the construction of the HCNetlas. Single-cell RNA sequencing data preannotated for cell type were used to build CGNs using the scHumanNet framework. HCNetlas is comprised of a comprehensive collection of these gene networks, representing various human tissues and cell types. (B) UMAP visualization of CGNs based on gene profiles, highlighting the major cell lineages, with node size representing the number of genes in each network. Major cell type “Other” in gray (major cell type abbreviations; B; B cells, Br; Brain cells, My; Myeloid cells, T; T cells). (C) UMAP plot displaying the interrelationship among the CGNs based on network gene profiles for major organs or tissue types. Each point represents a gene network associated with a specific organ or tissue type colored distinctly. The data underlying this figure can be found in https://doi.org/10.5281/zenodo.14522296. CGN, cell type-specific gene network; HCNetlas, human cell network atlas; UMAP, uniform manifold approximation and projection.

https://doi.org/10.1371/journal.pbio.3002702.g001

When constructing CGNs, the number of cells used can influence the efficacy of network inference. To investigate this aspect, we conducted an analysis of CGNs derived from various cell atlas data sets, specifically examining how the number of cells used for network inference correlates with the overall network size. Our findings indicated a clear trend: as the number of cells increases, there is a corresponding rise in both the node and edge counts within the inferred CGNs. However, this growth in network complexity tends to plateau once the cell count reaches approximately 1,000 (S1A Fig). The observed saturation point suggests that the inference of CGNs becomes substantially robust to the effects of sample size when the number of cells exceeds 1,000. Based on this insight, we focused our study on networks inferred from data sets comprising a minimum of 1,000 cells. This led to the generation of 198 CGNs, covering 25 tissues and including 61 distinct cell types (S1 Table). These networks form our newly established resource, HCNetlas, a catalog of human CGNs for healthy tissues.

To examine the interrelationships among the CGNs in our HCNetlas, we analyzed each CGN based on network gene profiles, subsequently visualizing these profiles in a reduced dimensional space. This analysis demonstrated a clear trend where CGNs corresponding to the same cell types exhibited a tendency to cluster together (Fig 1B), reinforcing the concept that these networks accurately capture and reflect the specificity inherent to each cell type. Notably, CGNs within the myeloid and B cell lineages showed remarkable coherence, in contrast to the T cell lineage CGNs, which exhibited greater heterogeneity. An interesting observation was the close proximity of innate lymphoid cell (ILC) and natural killer (NK) CGNs (S1B Fig), underscoring their lineage correlations [34]. However, CGNs related to the same tissue types generally did not demonstrate strong clustering with the exception brain tissue network nodes that displayed high similarity (Fig 1C), suggesting that cell type identity is a stronger determinant of network structure than tissue environment. This was further evidenced in the T cell lineage, including ILCs, NK cells, CD4⁺ T cells, and CD8⁺ T cells, where subsets exhibited coherence within cell types but not necessarily within tissue types (S1B and S1C Fig). This aligns with recent studies that emphasize tissue or sub-cell type-dependent variability in T cells [35–37]. These findings highlight the utility of HCNetlas as a potentially powerful tool for investigating cell type-specific gene functions.

To evaluate whether the reference CGNs of HCNetlas accurately reflect cell type-specific functions, we conducted tests using 2 distinct immune cell types from different lineages: B cells and T cells. The premise of this test was that if the HCNetlas CGNs are effective in representing cellular functions unique to each cell type, then genes for maintaining the identity and function of each cell type should demonstrate interconnectedness within their respective networks. As anticipated, our analysis showed that genes specifically annotated for either B cells or T cells by the GOBP [38] exhibited the highest within-group connectivity in their respective CGNs across various tissues (Fig 2A and 2B). This pattern of connectivity was further validated by comparing it with cell type marker genes as identified in the Azimuth database [13] and the CellMarker database [14] (S2A Fig). These findings underscore the ability of HCNetlas CGNs to capture and represent the unique functional characteristics inherent to specific cell types.

Fig 2. Cell type-specific functionality of HCNetlas CGNs.

(A, B) Bar graph illustrating the within-group connectivity of B cell-related (a) or T cell-related (b) GOBP genes in the respective CGNs. Connectivity is normalized by each CGN’s total node number. All tissues with over 0 value of normalized connectivity for both B and T cells are included. (C) Heatmap displaying the percentile rank of the top 15 hub genes in spleen CGNs, with values scaled per row, with color intensity indicating the expression level from low (blue) to high (red). Genes uniquely found in each cell type are highlighted in bold. (D, E) UpSet plots for 2 major CGNs, B cell CGN (D), and T cell CGN (E), representing the intersection of network genes across different tissues. The data underlying this figure can be found in https://doi.org/10.5281/zenodo.14522296.

https://doi.org/10.1371/journal.pbio.3002702.g002

Moreover, we investigated the network hub genes within each CGN, identified based on degree centrality (Fig 2C). For instance, in spleen CGNs, CD86, which is pivotal in B cell activation [39], emerged as top hub genes in the B cell CGN. Similarly, genes essential for T cell identity like CD2, CD4, and CD28 were among the top 15 hub genes in the T cell CGN [40]. Additionally, S100A8, S100A9, CD14, markers for monocytic myeloid-derived suppressor cells were prominent hubs in the monocyte CGNs [41]. These patterns of hub genes, significant due to their high degree centrality, were consistent across various tissues (S2B Fig), underlining the functional interpretability of these hub genes in the context of their respective cell types.

Lastly, to assess the tissue dependency of the HCNetlas CGNs, we compared CGN genes for each major cell type across different tissues. We observed limited overlap in CGN genes among different tissues within major cell types (Figs 2D, 2E, and S3), suggesting a convergence of networks across tissues within major cell lineages, aligning with findings from previous studies [10,35]. These observations indicate that while there are core gene networks characteristic of each cell type, tissue adaptation of CGNs is also evident, underscoring the complexity and diversity of cellular functions across different biological contexts.

The HCNetlas, with its collection of reference CGNs, presents a promising resource for dissecting the cellular specificity of disease genes. The majority of disease-associated genes identified to date have been derived from bulk tissue data, which often fails to specify the exact cell types involved in disease onset and progression. In this scenario, HCNetlas CGNs become instrumental in pinpointing the critical cell types at play. To ascertain the effectiveness of HCNetlas CGNs in disease-oriented research, we embarked on an investigation to determine if these CGNs capture and reflect the cell type specificity of various diseases. This involved conducting enrichment analyses on the CGNs using disease-associated genes sourced from 2 distinct databases: DisGeNET [16] and GWAS Catalog [17]. We began by ranking genes within each CGN based on network degree centrality, and then applied single-sample gene set enrichment analysis (ssGSEA) [18] and gene set variation analysis (GSVA) [19] to profile degree of association with each set of disease genes. Our analysis revealed a distinct pattern of congregation among CGNs corresponding to shared cell types, as determined by disease-association profiles (Figs 3A and S4A). This finding was particularly notable within cell types, whereas the convergence of networks corresponding to the same tissue types was less pronounced, indicating the specificity of cell types in the context of disease genetics.

We next evaluated the connectivity among genes associated with the same disease within CGNs across different tissues. Our hypothesis was that genes would exhibit more interconnectedness in the relevant cell types and tissues primarily responsible for diseases. This analysis aimed to elucidate the relationships between specific diseases and their associated cell types or tissues. While not all diseases we considered manifests cell type specificity, we noticed that CGNs predicted similar disease gene enrichment patterns in tissues such as the intestine and the liver (Fig 3A). A case in point is hepatitis-related terms, where genes associated with this condition showed the most significant within-group connectivity in liver CGNs of most major immune cell types (Fig 3B). Noteworthy was the pronounced within-group connectivity observed in both myeloid cell and T cell CGNs, highlighting the integral role of T cells in viral infectious diseases and the contribution of Kupffer cells (resident liver macrophages) to hepatitis [42]. This finding indicates that HCNetlas effectively identifies relevant cell types and tissues implicated in hepatitis. Furthermore, genes related to schizophrenia showed increased within-group connectivity across brain tissues, particularly in the primary motor cortex (M1) and MTG CGNs (S4B Fig). Additionally, AD-associated genes from the GWAS were found to be highly connected within the macrophage and myeloid CGNs (S4C Fig), consistent with the association between AD risk and macrophage transcriptional network [43]. These observations underscore the potential of HCNetlas CGNs as a valuable resource for uncovering intricate relationships between diseases and specific cell types or tissues, thereby enhancing our understanding of disease pathology at a cellular level.

Fig 3. Overview of cell type-resolved disease genetics using HCNetlas.

(A) Heatmap displaying the disease profiles of various cell types across different tissues, conducted with ssGSEA. Each column represents a CGN of HCNetlas, while each row corresponds to a disease gene set sourced from either DisGeNET or GWAS Catalog, ordered by hierarchical clustering. Color intensity indicates scaled ssGSEA enrichment values for the disease gene sets (rows), reflecting the degree of association of the CGN signature genes with each disease terms. (B) Bar graphs showing the within-group connectivity of genes associated with toxic hepatitis across different cell lineages, in various tissues. The bars are color-coded to represent different cell lineages; 15 terms and their combined genes were assessed from DisGeNET terms based on the key word search “hepatitis.” (C) Schematic representation and summary of the analytical framework used for comparing disease CGNs with reference CGNs from HCNetlas. The workflow illustrates the process of CGN inference from single-cell transcriptomes of disease samples and contrasts disease CGNs for specific immune cells against reference CGNs. The analysis includes (i) differential compactness, highlighting the difference in interconnectivity within disease-associated genes; (ii) differential hubness, showing the changes in hubness; and (iii) differential pathways, contrasting pathway associations between disease and healthy states based on enrichment for CGN signature genes. The data underlying this figure can be found in https://doi.org/10.5281/zenodo.14522296. CGN, cell type-specific gene network; HCNetlas, human cell network atlas; ssGSEA, single-sample gene set enrichment analysis.

https://doi.org/10.1371/journal.pbio.3002702.g003

HCNetlas, having proven its functional and biological relevance, is posited to be an effective reference for network analyses in disease studies. To enhance their utility, we have developed a suite of network analysis methodologies (Fig 3C) and applied them to investigate various diseases, showcasing the adaptability of HCNetlas CGNs.

Given that the majority of the CGNs provided by HCNetlas are derived from immune cells, this resource would be particularly valuable for studying immune disorders such as autoimmune diseases. To evaluate capability of HCNetlas to identify the specific immune cell types where disease genes have an impact on pathogenesis, we focused our research on SLE which is a chronic autoimmune disorder characterized by elusive pathogenesis, genetic susceptibility, and clinical heterogeneity [44]. For constructing disease CGNs for SLE, we manually annotated scRNA-seq data from 38 SLE patients [20] using canonical markers (Fig 4A). These disease CGNs were then compared with the blood cell CGNs from HCNetlas, providing insights into the cell type specificity underlying SLE pathogenesis.

Fig 4. Cell type-resolved disease genetics for SLE.

(A) UMAP plot representing the interrelationship among immune cells. (B) Bar chart showing normalized interconnectivity among SLE-susceptible genes in both reference and disease CGNs for various cell types. (C) ROC curves for retrieval of ISGs by network hubness in disease myeloid CGN and reference myeloid CGN. (D) Comparison of AUROC values with CGNs for various cell types, contrasting the reference and disease CGNs to assess prediction capability. (E) Left panel: Distribution of expression levels for 131 gain-of-hubness genes in myeloid cells. Right panels: Correlation analysis of expression levels of the 131 gain-of-hubness genes with the SLEDAI. (F) Same as for (E) except using up-regulated DEGs. The data underlying this figure can be found in https://doi.org/10.5281/zenodo.14522296. AUROC, area under the ROC curve; CGN, cell type-specific gene network; DEG, differentially expressed gene; ISG, interferon-stimulated gene; ROC, receiver operating characteristic; SLE, systemic lupus erythematosus; SLEDAI, Systemic Lupus Erythematosus Disease Activity Index; UMAP, uniform manifold approximation and projection.

https://doi.org/10.1371/journal.pbio.3002702.g004

Leveraging the principle that increased network compactness among disease-associated genes within a CGN indicates their significant role in the pathogenesis for that cell type, we assessed the involvement of major immune cell types in SLE. We applied a set of SLE-susceptible genes (S2 Table), gathered from the KEGG pathway database, to both disease CGNs and reference CGNs. Our analysis revealed that network compactness in myeloid cells and B cells is significantly greater in the disease CGN compared to the reference CGN (Fig 4B). This suggests that SLE-susceptible genes are critically involved in the disease progression primarily through myeloid cells and B cells.

Considering that genes associated with SLE predominantly exert their effects through myeloid cells, we prioritized genes for SLE based on network centrality within both disease and reference CGNs specifically pertaining to myeloid cells. Aligning with previous studies that emphasize the increased expression of type 1 interferon (IFN) and ISGs in SLE patients [20,45,46], we evaluate the prediction of SLE-associated genes based on retrieval rate of ISGs (S3 Table) using the ROC curve. Consistent with the greater network compactness of SLE-susceptible genes in the disease myeloid CGN relative to the reference myeloid CGN, our results showed a significantly improved prediction of ISGs in the disease myeloid CGN when compared to the reference myeloid CGN (Fig 4C). Likewise, for other cell types, disease CGNs demonstrated improved predictions of ISGs compared to the reference CGNs (Fig 4D and S4 Table). Taken together, these findings underscore the critical role of myeloid cells in the initiation and progression of SLE, corroborating previous research that highlights the link between SLE and myeloid cells [47–49].

Next, we hypothesized that gain-of-hubness genes for myeloid cells could effectively differentiate diseased myeloid cells from their healthy counterparts. To test this hypothesis, we initially identified a set of 131 gain-of-hubness genes with statistical significance (S5A Table). We then examined the distribution of their expression level between disease-state myeloid cells and their corresponding healthy controls. Our observations revealed a significant disparity between these 2 distributions, affirming the potential of these 131 gain-of-hubness genes to distinguish diseased myeloid cells (Fig 4E, left panel). Furthermore, we observed a positive correlation between the expression levels of these gain-of-hubness genes and the Systemic Lupus Erythematosus Disease Activity Index (SLEDAI) scores, albeit with limited statistical power due to the small sample size (Fig 4E, right panel). This correlation indicates that the expression patterns of these 131 gain-of-hubness genes are not only distinctive of diseased states but may also reflect the severity of SLE in patients. In contrast, the up-regulated DEGs in disease-state myeloid cells (S5B Table) did not demonstrate the capability to either differentiate diseased myeloid cells (Fig 4F, left panel) or correlate with SLEDAI scores (Fig 4F, right panel). These outcomes imply that collections of CGNs for healthy tissues, such as those provided by HCNetlas, are apt references for identifying disease states.

HCNetlas offers an extensive collection of CGNs for brain tissue, making it a valuable resource for investigating neurological disorders. Alzheimer’s disease (AD), a widespread neurodegenerative condition known for its progressive impact on behavior and cognitive functions, is one such area where HCNetlas can be particularly useful. To study AD more closely, we have developed disease CGNs using single-nucleus RNA sequencing (snRNA-seq) data from the prefrontal cortex of AD patients [25]. These disease CGNs were compared with HCNetlas CGNs, which were derived from the primary motor cortex (M1). This comparison enables a detailed analysis of alterations in the gene network that is associated with AD, facilitating a deeper understanding of the disease progression and its impact on brain function.

To identify the primary cell types impacted by AD-associated genes compiled from various sources (Methods, S6 Table), we employed differential compactness analysis. This analysis revealed that AD-associated genes exhibit a high degree of interconnectivity within the reference CGNs for both inhibitory and excitatory neurons (Fig 5A), suggesting that these genes predominantly function in these neuron types. Interestingly, we observed that the disease CGNs for inhibitory and excitatory neurons displayed notably lower network compactness scores compared to their reference counterparts (i.e., loss-of-compactness). This significant decrease in network compactness within the diseased neurons points to a loss of connections among AD-associated genes. Such a loss in the diseased state of inhibitory and excitatory neurons could be a critical factor in the pathogenesis of AD, indicating a disruption in the intricate gene networks that underlie normal neuronal function.

Fig 5. Cell type-resolved disease genetics for AD.

(A) Bar graph depicting the normalized edge count among AD-associated genes in reference and disease CGNs for various neurological cell types. (B) Bar graph showing the number of validated AD genes predicted by differential expression or differential hubness across the 4 neurological cell types. Statistical significance of overlap is shown for each gene sets (*P P P C–F) Ten most differentially associated KEGG pathways with CGN signature genes: Differentially associated pathways for inhibitory neuron CGN signature (C), inhibitory neuron expanded CGN signature (D), excitatory neuron CGN signature (E), and excitatory neuron expanded CGN signature (F). The data underlying this figure can be found in https://doi.org/10.5281/zenodo.14522296. AD, Alzheimer’s disease; CGN, cell type-specific gene network.

https://doi.org/10.1371/journal.pbio.3002702.g005

We then focused on prioritizing genes for AD by either differential expression or differential hubness between healthy and diseased states across each cell type (S7 Table). In alignment with the identified cell type specificity for AD, both inhibitory and excitatory neurons demonstrated a more accurate prediction of AD-related genes when analyzed for differential hubness rather than differential expression (Fig 5B, hypergeometric test P-value S5A Fig), indicating that these 2 approaches are complementary to each other in identifying key genes associated with AD.

To delve into the molecular mechanisms implicated in AD pathogenesis within inhibitory and excitatory neurons, we carried out a differential pathway analysis. This analysis was based on CGN signatures of these neurons, focusing on the top 10 genes ranked by hubness (S5B Fig). As anticipated, our analysis revealed that pathways associated with AD and other related neurodegenerative diseases, such as Parkinson’s disease and Huntington’s disease, were among those most prominently exhibiting a reduced association, or “loss-of-pathway,” in inhibitory neurons (Fig 5C). In addition to these, we identified several other pathways that exhibited loss-of-pathway in inhibitory neurons, and these findings were validated through a literature survey. The pathways that were validated to be associated with AD included oxidative phosphorylation [50], thermogenesis [51], non-alcoholic fatty liver disease [52,53], diabetic cardiomyopathy [54,55], amyotrophic lateral sclerosis [56], and prion disease [57]. Significantly, our analysis identified that the pathway related to the cholinergic synapse was the most notably increased in diseased inhibitory neurons. This finding is also relevant given the known involvement of the cholinergic signaling in AD [58]. We also performed our analysis using an expanded CGN signature that includes their network neighbors, which confirmed the initial list of top loss-of-pathways (Fig 5D). This reaffirms the importance of these pathways in the pathogenesis of AD, highlighting their potential roles in the disease’s progression and impact within inhibitory neurons.

In our differential pathway analysis using CGN signatures for excitatory neurons, we predominantly observed gain-of-pathways. These are pathways showing increased activity in AD compared to the healthy state, findings which are substantiated by literature evidence (Fig 5E). For instance, the ErbB signaling pathway is known to mediate amyloid-β (Aβ)-induced neurotoxicity [59], and HIF-1 (hypoxia-inducible factor-1) signaling has been found to increase Aβ generation [60]. Additionally, a similar analysis with expanded CGN signatures for excitatory neurons revealed loss-of-pathways akin to those identified in inhibitory neurons (Fig 5F). Among these findings, the MAPK signaling pathway emerged as the most prominent gain-of-pathway. This is in alignment with previous research demonstrating that MKP-1, a crucial negative regulator of MAPKs, can reduce Aβ generation and alleviate cognitive impairments in AD models [61], thereby validating our observation.

The tumor immune microenvironment has become increasingly recognized as a key hallmark of cancer. Considering this, we hypothesized that HCNetlas CGNs for immune cells could be instrumental in identifying cancer-associated genes that primarily function within the immune microenvironment. Focusing on lung cancer, we constructed disease CGNs for major immune cell types using single-cell transcriptome data derived from tumor tissues of lung cancer patients [28]. Through differential hubness analysis, compared to reference CGNs of corresponding cell types, we pinpointed gain-of-hubness genes predominantly in T cells and myeloid cells, many of which are known to be associated with lung cancer (Fig 6A and S8A Table). Interestingly, only a few gain-of-hubness genes were common across multiple immune cell types, suggesting a specific functional role of cancer-associated genes in T cells and myeloid cells. This analysis also revealed that differential hubness was more effective in identifying lung cancer-associated genes than the traditional differentially expressed genes (DEGs) analysis (Fig 6B and S8B Table). Notably, many up-regulated DEGs shared among all immune cell types included very few validated lung cancer genes. When assessing loss-of-hubness genes, a similar trend was observed: fewer candidates but with more specificity to cell types compared to down-regulated expression in disease (Fig 6C and 6D). T cell-specific loss-of-hubness particularly retrieved a significant number of known lung cancer genes. Additionally, we found supportive literature evidence for the proposed cell type of action for these validated cancer-associated genes identified through differential hubness analysis (S8 Table). This suggests that HCNetlas is effective in predicting genes associated with cancer specifically within immune cell types.

Fig 6. Identifying genes that contribute to lung cancer through immune cells using HCNetlas.

(A–D) UpSet plots for predicted genes for lung cancer by gain-of-hubness (A), up-regulated DEGs (B), loss-of-hubness (C), or down-regulated DEGs (D), representing the intersection of predictions across different cell types. Orange bar indicates the number of lung cancer genes validated by various databases such as CancerMine, IntOGen, and cancer gene consensus. (E) Bar graph showing the number of ICMs retrieved by both gain-of-hubness genes and up-regulated DEGs across different immune cell types. Statistical significance of overlap is shown for each gene sets (*P P P F, G) Kaplan–Meier survival curve analysis for cancer patients from TCGA cohort (TCGA-LUSC and TCGA-LUAD), stratified by the enrichment score of gain-of-hubness genes (F) or up-regulated DEGs (G). The graph shows the overall survival probability over time, with patients categorized into high and low quantiles based on the GSVA score. Significance of survival rate difference between the upper and lower quantile expression groups were evaluated using the log-rank test. The data underlying this figure can be found in https://doi.org/10.5281/zenodo.14522296. DEG, differentially expressed gene; HCNetlas, human cell network atlas; GSVA, gene set variation analysis; ICM, immune checkpoint molecule.

https://doi.org/10.1371/journal.pbio.3002702.g006

Further evaluation focused on immune checkpoint molecules (ICMs), which are pivotal in antitumor immunity [62,63]. We anticipated an increase in network centrality and expression of ICMs in tumor-derived immune cells. Confirming our hypothesis, genes identified through differential hubness analysis were more effective in detecting ICMs, particularly within T cells and myeloid cells, compared to differential expression analysis (Fig 6E). This finding underscores the advantage of network-based analyses in pinpointing crucial genes in cancer immunology. Additionally, using The Cancer Genome Atlas (TCGA) lung cancer data, we explored the prognostic value of these genes. We found that the gene expression profile association score, calculated using GSVA [19], for the set of gain-of-hubness genes in each tumor sample was predictive of clinical outcomes (Fig 6F), unlike up-regulated DEGs (Fig 6G).