Scalable cell-specific coexpression networks for granular regulatory pattern discovery with NeighbourNet

Coexpression network

NNet estimates coexpression at the level of individual cells by combining dimensionality reduction with local regression. To measure cell-specific gene coexpression, NNet first uses principal component analysis (PCA) to embed gene expression data into a low-dimensional space (Fig. 1A). Within this space, each cell finds its k-nearest neighbors (KNNs), and a regression model is fitted to predict the expression of a response gene using the PCs. Within the neighborhood of cell n, the coexpression level between the response gene p and a gene q contributing to PC computation (a predictor gene), denoted CSN_npq, is defined as how strongly the gene q influences the prediction of p through PCs given the fitted model.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

NNet workflow for inferring cell-specific coexpression networks and integrating prior knowledge. (A) Preprocessing: Single-cell gene expression data are subject to PCA to capture the major variation. A low-rank approximation (LRA) then reconstructs the expression matrix using the PCs. A weighted k-nearest neighbor (KNN) graph that defines each cell’s local neighborhood in the PC space is constructed. (B) Neighborhood regression: For each cell, NNet performs regression within the cell’s KNN using the LRA-derived gene expression as the response and PCs as predictors. This process quantifies coexpression between genes by measuring how “predictor genes” (those used to embed PCs) contribute to predict each response in the PC space. Repeating this for multiple response genes yields cell-specific coexpression networks (CSNs). (C) Network ensemble construction: Collecting all CSNs form a network ensemble, which is a cell × response × predictor 3D array that can be diced and sliced. (D) Downstream analysis: (D1) We can slice the network ensemble to examine CSNs, analyzing cellular variation and gene regulation dynamics through coexpression. (D2) Nonnegative matrix factorization (NMF) identifies soft clusters of cells with similar network structures. Aggregating CSNs by these clusters yields meta-networks that summarize shared coexpression patterns. Same analysis strategy can be applied to the gene dimension corresponding to TFs, producing TF modules that coregulate common programs. Aggregating TFs by module yields simplified meta-TF–target networks with improved interpretability. (D3) We adapted the NicheNet framework to integrate gene regulation and signaling interaction databases from OmniPath, constructing integrated prior knowledge networks (PKNs). Annotating CSNs with these PKNs transforms them into contextualized gene regulatory networks (GRNs). In addition, PKNs enable the inference of upstream signaling pathways (USPs) for each CSN by tracing signal transduction paths (receptor–TF–target) inferred based on the contextualized TF–target interactions.

By repeating this regression procedure across all response genes, NNet constructs a weighted coexpression network CSN_n.. for every individual cell, represented by a matrix of coexpression values between predictor and response genes for that cell (Figs. 1B,C). As NNet is regression-based, switching the roles of response and predictor genes yields different coexpression values for the same gene pair. The choice of responses depends on context, with the general rule of keeping the number of responses smaller than the number of predictors to reduce computational cost.

Embedding network ensemble

The collection of CSNs builds a network ensemble, which is a cell × response × predictor three-dimensional array (a “cube”) where each slice along the cell dimension corresponds to a CSN. While this per-cell resolution is powerful, it can also be helpful to extract higher-level structures that summarize thousands of such CSNs. To achieve this, NNet employs nonnegative matrix factorization (NMF) to identify coherent coexpression patterns across the ensemble, grouping cells that share similar network structures into soft clusters and summarizing them as corresponding meta-networks (Fig. 1D1,2).

Each soft cluster of cells is represented by a matrix H, where each column is a weighting vector defining cluster i, and each entry indicates the degree of membership of cell n to that cluster. NMF learns these soft clusters based on cell–cell similarity in their networks, and the corresponding meta-network for cluster i is computed as a weighted average of the individual CSNs: Our NMF approach extracts multiple soft clusters and meta-networks in sequential order, such that meta-network 1 captures the most significant coexpression pattern present in the data set.

The same strategy can be applied to learn soft clusters of genes, summarizing slices of the network ensemble along either the response or predictor dimension. When applied to TFs, each soft cluster represents a TF module, which is a group of TFs that act together to regulate a common transcriptional program. Aggregating TFs according to their learned module yields meta-TF–target CSNs with only a few meta-TFs.

Signaling inference

In parallel, we incorporate prior knowledge by mapping observed coexpression edges to curated regulatory and signaling interaction databases. This step produces context-specific gene regulatory networks (GRNs) with regulatory directories and allows for the inference of upstream signaling pathways (USPs) that connect ligand–receptor interactions to transcriptional responses through contextualized TF–target coexpression (Fig. 1D3). The output includes interpretable networks, tools for visual exploration of networks’ regulatory structures, and quantification of active receptor signaling across individual cells.

Data and setting

We analyzed two data sets. The first was the PBMC 3K data set from 10x Genomics, often used as a benchmark. The second data set was a Perturb-seq experiment (Papalexi et al. 2021), in which the THP-1 monocyte cell line was subjected to pooled CRISPR screening to characterize the molecular networks regulating PD-L1 expression (Table 1).

For the PBMC data set, our aim was to perform a sanity check to assess whether TF activity inference based on NNet-derived coexpression is biologically meaningful. Specifically, we analyzed 288 TFs with sufficient target coverage according to the CollecTRI gene regulation database (Müller-Dott et al. 2023) (Methods). Using NNet, we computed each TF’s coexpression with 4116 potential target genes (i.e., any gene with a known TF regulator in CollecTRI). From these coexpression profiles, we derived CoETAS and compared them with ETAS generated directly from target gene expression. We then evaluated which scores better preserved cell-type clustering and better highlighted known cell-type-specific TFs. This comparison serves as a validation step, showing that NNet’s coexpression measure provides sufficient biological signal to capture meaningful regulatory variation.

For the Perturb-seq data set, which includes perturbations on 5 TFs across different subsets of cells, we again used NNet to measure coexpression between perturbed TFs and 5151 target genes. Subsequently, we calculated ETAS and CoETAS for each of the 5 perturbed TFs across all cells. Because the data set provides precise information on TF perturbations in each cell, it allowed us to directly assess how well the inferred activity scores distinguish perturbed from unperturbed cells. This, in turn, provides a hint of whether NNet’s coexpression measure carries accurate information about TF regulatory status.

NNet coexpression-based TF activity scores better preserves cell clusters in PBMC data

We annotated the data set by first performing cell clustering based on the expression of the 2000 most variable genes, followed by cell type assignment using cluster-specific marker genes (Fig. 2B1). Using NNet, we constructed CSNs, which enabled novel analyses of cellular variation based on network topology, such as connectivity (i.e., node degree). A UMAP embedding of cells based on the connectivity profiles of 288 TFs within these networks largely preserved the separation of cell types (Fig. 2B2).

Using both target gene expression and NNet TF–target coexpression, we computed TF activity scores and assessed how well these scores preserve cell type clusters. UMAP embeddings showed that CoETAS (Fig. 2B3) better captured known cell types compared to ETAS (Fig. 2B4, Supplemental Fig. S1A). Quantitative evaluations using the adjusted Rand index and the silhouette index further confirmed that CoETAS produced higher-quality clustering (Fig. 2C). As a control, clustering based on TF expression performed worse than TF activity. Additional evaluations using clusters derived from low-rank approximations of TF expression (Supplemental Fig. S1B) ruled out the possibility that the advantage of CoETAS stemmed from the PCA regression used to measure coexpression. These results highlight the ability of NNet coexpression in capturing cell-specific regulatory activity.

NNet coexpression profiles detect more meaningful cell-type-specific gene regulation in PBMC data

We then assessed the biological relevance of differentially activated TFs (DATs) identified by activity scores in different cell type clusters. We focused on B cells, which formed the most distinct cluster within the PBMC data (Fig. 2B). We conducted t-tests on activity scores to extract DATs in B cells. The t-tests on ETAS and CoETAS revealed largely different sets of DATs (Fig. 2D). Whereas most DATs identified by both methods were differentially expressed by B cells, several DATs identified solely by ETAS were actually downregulated or not expressed at all by the B cell population. A key example was TBX21, a master regulator of T and NK cells (Szabo et al. 2000; Knox et al. 2019), which was deemed active in B cells by ETAS. In contrast, CoETAS correctly highlighted the relevance of TBX21 in T cells and NK cells (Fig. 2E). Thus, relying solely on target expression for TF activity inference can yield spurious results, whereas NNet coexpression provides more accurate evidence for genuine regulatory activity.

NNet coexpression-based TF activity scores distinguish perturbation effects more effectively in the Perturb-seq data

We further assessed the robustness of CoETAS by analyzing TF perturbation (knockdown) experiments in the THP-1 cell line. A UMAP embedding of cells showed variable responses to five TF perturbations, with STAT1 knockdown producing the strongest effect, as perturbed cells were clearly separated from the control cluster (Fig. 2F). For each perturbed TF, we computed CoETAS and ETAS using decoupleR, and used area under the curve (AUC) to evaluate their ability to distinguish perturbed cells (expected lower activity; negative class) from controls (expected higher activity; positive class). On average, CoETAS outperformed both ETAS and raw TF expression, achieving higher AUC scores (Fig. 2G). The largest performance gain was observed for SMAD4, where CoETAS clearly identified reduced activity in SMAD4-perturbed cells (Fig. 2H), whereas ETAS patterns were less distinct (Supplemental Fig. S1C1). In contrast, CoETAS performed worst for STAT2, which is likely due to the weak perturbation effect suggested by the close proximity of STAT2-perturbed cells to controls in the UMAP. Nevertheless, CoETAS still highlighted a subcluster enriched for STAT2-perturbed cells (Supplemental Fig. S1C2).

We performed a similar analysis on an independent Perturb-seq data set (Dixit et al. 2016) containing more TF perturbations (Supplemental Results; Supplemental Fig. S1) and confirmed the advantage of CoETAS over ETAS in detecting perturbation effects.

Overall, our findings show that incorporating NNet coexpression significantly enhances the accuracy and biological relevance of TF activity inference, underscoring its potential in effectively capturing gene regulation activities. Additional analysis comparing NNet with two other cell-specific methods, oCSN (Dai et al. 2019) and LocCSN (Wang et al. 2021), revealed that NNet substantially improves CSN construction by better preserving cell-type identity (Supplemental Results; Supplemental Fig. S7).

Data and setting

We used the lineage-negative (Lin⁻) bone marrow cell scRNA-seq data from (Pellin et al. 2019), which was designed for profiling the early human hematopoietic landscape (Table 1). With TF and target gene sets acquired from our integrated prior knowledge network (PKN) of gene regulation (Methods), NNet analysis was performed to measure coexpression between 805 TFs and 4600 potential targets.

In the following sections, we perform meta-network analysis and meta-TF analysis separately to study two distinct major lineages: the erythroid lineage and the mononuclear phagocyte (MP) lineage. Meta-network analysis is used to examine transcriptional program alterations during erythropoiesis, whereas meta-TF analysis is applied to identify the key TF module that potentially coordinates MP differentiation.

NNet meta-network analysis reveals abrupt rewiring of coexpression network during erythrocyte maturation

We identified 20 meta-networks, each characterized by a set of cell weights indicating the cell populations (i.e., soft cluster) they represent (Supplemental Fig. S2A; see Methods for a detailed definition of cell weights). These meta-networks captured key transcriptional shifts during hematopoiesis. Specifically, the first four meta-networks delineated the primary myelopoietic lineages towards erythroid (Fig. 3B) and MP fates (Supplemental Fig. S2B,C). For example, meta-network 1, primarily composed of networks of megakaryocyte–erythroid progenitors (MEPs), was associated with early lineage commitment characterized by HDAC7-driven suppression of MYC that releases erythroid differentiation blockage (Fig. 3C) (Delgado and León 2010; Jayapal et al. 2010; Wang et al. 2020). In contrast, meta-network 2, which exhibited increased involvement during terminal erythropoiesis, highlighted a shift towards network modules related to hemoglobin synthesis (e.g., KLF1-mediated ABCB10 expression) (Tallack et al. 2010) and erythrocyte membrane protein synthesis (e.g., TAL1-mediated RHD and ERMAP expression) (Kassouf et al. 2010). Furthermore, the cell weights on meta-network 2 sharply localized to erythroid progenitors (ErPs) from MEPs, suggesting that erythrocyte maturation may involve abrupt, temporally distinct transcriptional reprogramming. Our results recapitulate a well-established TF regulation mechanism known as the “GATA switch,” in which the repression of GATA2 facilitates a surge in GATA1 expression: a crucial step in erythrocyte maturation (Suzuki et al. 2013; Bresnick et al. 2018).

NNet network topology of individual cells better characterizes temporal dynamics of GATA regulation compared to gene expression

Focusing on the GATA genes (GATA1 and GATA2), we explore how NNet network topology analysis provides an enhanced characterization of dynamics of transcriptional program compared to conventional gene expression analysis. We calculated the connectivity of GATA genes within CSNs and plotted this connectivity against the diffusion pseudotime of cells (Fig. 3D1) (Haghverdi et al. 2016). This analysis showed that GATA2 reached its peak connectivity at the MEP stage, which then sharply declined following the transition to the ErP stage. As GATA2’s connectivity dropped to zero, GATA1 began to show a significant increase in connectivity throughout maturation, emphasizing a critical temporal coordination of GATA2 and GATA1 during erythropoiesis. In contrast, plotting the expression of GATA genes against pseudotime lacks critical details about the dynamics of GATA regulation, such as the peak of activity and the switching point between transcriptional states (Fig. 3D2).

NNet meta-TF analysis identifies key transcriptional programs that potentially direct CD141⁺ dendritic cell differentiation

We proceed to illustrate the power of meta-TF analysis, focusing on the identification of key meta-TFs and their regulons associated with MP differentiation, a process distinct from erythropoiesis. Each meta-TF represents a TF module (i.e., a soft cluster) defined by a set of TF weights. Although analyzing a single TF’s coexpression often provides only a partial view of its function, interpreting TFs within the context of their modules offers an integrative view of their role as a regulatory program.

In each CSN, TFs are aggregated according to modules to analyze meta-TF level coexpression with targets (Supplemental Fig. S2D). By examining the connectivity of meta-TFs within each cell, we identified meta-TF 7 as the first one that showed strong connectivity across the MP lineage, particularly in granulocyte–monocyte progenitors (GMPs) (Fig. 3E; Supplemental Fig. S2E). This suggests a potential role for meta-TF 7 in MP differentiation. To functionally characterize meta-TF 7, we performed a Gene Ontology (GO) over-representation analysis on its module composition. This analysis revealed that meta-TF 7 was closely associated with immune cell differentiation (Fig. 3F), particularly the differentiation of dendritic cells (DCs), which belong to the MP lineage. Notably, the DC differentiation term was enriched by BATF3, a key TF that is essential for the development of conventional type 1 DCs (cDC1) (Hildner et al. 2008; Satpathy et al. 2012). This clinically significant DC subtype plays a critical role in anti-cancer immunity (Böttcher et al. 2018), with extensive research dedicated to developing protocols for its in vitro differentiation (Elahi et al. 2022; Rosa et al. 2022). We then sought to identify and comprehend potential downstream targets of meta-TF 7 to gain a deeper understanding of its function. Figure 3G illustrates the changes in target gene coexpression with meta-TF 7 during MP differentiation. By clustering targets based on these coexpression values, we identified two major regulons that exhibited staggered waves of activity (Supplemental Fig. S2F).

The first wave of activity, represented by regulon 1, included IRF8, KLF4, and CSF1R, which are critical markers indicative of MP fate commitment (Kurotaki et al. 2014; Tussiwand and Gautier 2015; Combes et al. 2021). The second wave, represented by regulon 4, featured BATF3 alongside other key markers of cDC1 identity, including ID2, a TF that delineates the cDC1 lineage (Jackson et al. 2011); THBD, which encodes the CD141 surface marker (Anastasiou et al. 2012); and CADM1, which is expressed exclusively by cDC1 (Collin and Bigley 2018). These findings further support the role of meta-TF 7 in coordinating the transcriptional programs driving cDC1 differentiation. Notably, our analysis of marker expression revealed that cDC1 constituted only a small population at the tip of the MP lineage, making it challenging to identify through clustering analysis alone. This underscores the effectiveness of our meta-TF analysis in disentangling relevant TF-target interactions within rare cell populations.

In conclusion, our analysis of the early hematopoiesis atlas demonstrates NNet’s capability to analyze cellular variation through coexpression, particularly providing insights into the dynamics of gene regulation throughout cell differentiation. By meta-network (meta-TF) analysis, we successfully identified principal coexpression patterns, TF modules and regulons at the individual cell level without the need for cell clustering.

Data and setting

We acquired the SCLC atlas from Chan et al. (2021), which includes scRNA-seq data of both tumor and tumor-adjacent normal tissues (Supplemental Fig. S3A; Table 1). Focusing on the macrophage subset (Fig. 4B), we conducted a differential expression (DE) analysis comparing macrophages from tumor and normal samples. Among the top 50 genes upregulated by tumor macrophages, we identified 29 genes that are targets of known TFs according to our PKN of gene regulation (Fig. 4C). The coexpression of these 29 genes with 900 TFs obtained from the PKN was then measured using NNet. We applied meta-network analysis to recover the major coexpression patterns and subsequently utilized USP inference to investigate how signals are transmitted to the 29 genes through TFs across different macrophage landscapes.

NNet annotated coexpression networks capture key signaling pathways defining macrophage interactions within TME

We focused on the top three meta-networks and projected the soft clusters they represent onto a UMAP embedding of macrophages (Fig. 4D). These soft clusters distinctly localized to different macrophage populations, suggesting that each population is governed by a unique transcriptional program. To further explore signaling differences among these populations, we visualized the annotated coexpression networks of the three most representative cells (those with the highest cell weight on each meta-network), highlighting a spectrum of macrophage identities (Supplemental Fig. S3B).

On one end of the spectrum, cell α, representing meta-network 1, showed a high level of connection across diverse regulatory modules (Fig. 4E). The USP analysis revealed an enrichment of proinflammatory signaling in cell α, with NFKB1 acting as the central TF that mediated signals from receptors like IL1R (IL1R1), TNFR1 (TNFRSF1A), and TLR4 (TLR4). This observation suggests the activation of the well-known NF-kB axis, which strongly supports the M1 proinflammatory identity (as opposed to the M2 anti-inflammatory identity) of cell α (Hagemann et al. 2009; Yu et al. 2020).

On the other end of the spectrum, cell γ, representing meta-network 3, displayed a complete loss of proinflammatory signaling (Supplemental Fig. S3C1). Instead, it established ENO1-mediated signaling to stimulate SPP1 expression, with SRB1 (SCARB1) as the primary receptor and SRC likely acting as the transducer relaying the signal to ENO1. SPP1⁺ macrophages represent a distinct subtype of TAM that has been reported to interact with fibroblasts and anti-inflammatory CD8⁺ T cells, driving tumor growth, metastasis, and immunosuppression (Qi et al. 2022; Bill et al. 2023). While the roles of SRB1, SRC, and ENO1 in TAM have been independently reported, their interaction in promoting SPP1 expression presents a potential regulatory mechanism that warrants further investigation (Dwyer et al. 2017; Plebanek et al. 2018; Liang et al. 2023).

In between cell α and γ, cell β, representing meta-network 2, exhibited an intermediate transcriptional state with an increase in the activity of the IL4R receptor, which is a critical switch for M2 state activation (Supplemental Fig. S3C2). We hypothesise that cell β and the macrophage population associated with meta-network 2 exist in a highly plastic state, rendering them susceptible to reprogramming and valuable targets for therapeutic treatment. To better understand this transitional state, our next analysis investigates the key signaling interactions and their consequences among these transient macrophages in the TME.

NNet receptor activity analysis highlights activation of RBPJ-mediated NOTCH signaling during transitional macrophage state

NNet’s USP inference involves calculating each receptor’s potential activity across target genes for each cell, summarized as activity scores. Leveraging this inference, we explored the receptors that exhibit differential signaling in transient macrophages, represented by cell β, compared to M1 macrophages (e.g., cell α) and the TAM population (e.g., cell γ). Specifically, we summarized receptor activity within each macrophage population by a weighted sum (using importance scores) of cell-specific receptor activities (Methods). Our comparison revealed that NOTCH (NOTCH2/3/4) receptors were exclusively active in the transient macrophage population associated with meta-network 2 (Fig. 4F). Endothelial cells, expressing the NOTCH ligands DLL and JAG, emerged as the principal source of NOTCH signaling (Supplemental Fig. S3D). This finding aligns with previous research showing that NOTCH signaling can influence macrophage polarization towards the M1 or M2 phenotype, although its precise direction remains a debated topic (Xu et al. 2012; Foldi et al. 2016; Palaga et al. 2018; Chen et al. 2023).

To clarify the role of NOTCH in our lung cancer setting, we examined TF connectivity in relation to NOTCH2 activity (Fig. 4G). RBPJ emerged as the top TF whose connectivity highly correlated with NOTCH2 activity (Fig. 4H), which is consistent with its established role as a key mediator of NOTCH signaling and a frequent target for blocking the pathway (Friedrich et al. 2022). Our findings partially recapitulate those of Franklin et al. (2014), who reported that RBPJ-mediated NOTCH signaling is crucial for the final differentiation of TAMs from tumor-infiltrating monocytes. However, it remains unclear whether this pathway also governs the macrophage transition from the M1 to TAM phenotype. To address this question, we next investigated the downstream targets of RBPJ to characterize the functional consequences of NOTCH signaling.

RBPJ coexpression pinpoints downstream targets of NOTCH signaling and its role in converting macrophages to TAM

We performed a separate NNet analysis to measure RBPJ’s coexpression with 5011 target genes. By calculating the correlation between NOTCH2 activity (as shown in Fig. 4G) and each target’s coexpression with RBPJ across individual cells, we identified a set of genes that were highly associated with NOTCH-RBPJ signaling (Fig. 4I). According to prior knowledge, key genes triggered by this signaling include MAF and DAB2, which are important regulators of M2 polarization (Liu et al. 2020; Marigo et al. 2020), as well as the cell surface receptors MRC1 (also known as CD206) and CD209, which are defining markers for M2-like TAMs (Wang et al. 2024). Whereas Franklin et al. (2014) suggested that RBPJ-dependent TAMs are MRC1⁻, our findings indicate that RBPJ may induce MRC1⁺ TAMs at a transient state. Functionally, we also noted that NOTCH-RBPJ signaling triggered macrophages to secrete chemoattractants such as CCL13, which may recruit immune cells to the TME, thereby modulating the immune landscape (Supplemental Fig. S3D).

Taken together, our final case study showcases NNet’s ability to capture signaling interactions that precede gene regulation in individual cells. Through this approach, we identify that NOTCH-RBPJ signaling becomes active at an intermediate macrophage state, positioned between the M1 proinflammatory and M2-like TAM phenotypes. In this transitional phase, NOTCH-RBPJ may modulate re-education signals from the TME, facilitating the conversion of M1 macrophages into TAM.