Single-cell biological network inference using a heterogeneous graph transformer

Single-cell multi-omics (scMulti-omics) allows the quantification of multiple modalities simultaneously to capture the intricacy of complex molecular mechanisms and cellular heterogeneity. Existing tools cannot effectively infer the active biological networks in diverse cell types and the response of these networks to external stimuli. Here we present DeepMAPS for biological network inference from scMulti-omics. It models scMulti-omics in a heterogeneous graph and learns relations among cells and genes within both local and global contexts in a robust manner using a multi-head graph transformer. Benchmarking results indicate DeepMAPS performs better than existing tools in cell clustering and biological network construction. It also showcases competitive capability in deriving cell-type-specific biological networks in lung tumor leukocyte CITE-seq data and matched diffuse small lymphocytic lymphoma scRNA-seq and scATAC-seq data. In addition, we deploy a DeepMAPS webserver equipped with multiple functionalities and visualizations to improve the usability and reproducibility of scMulti-omics data analysis.

Introduction

Single-cell sequencing, such as single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq), reshapes the investigation of cellular heterogeneity and yields insights in neuroscience, cancer biology, immuno-oncology, and therapeutic responsiveness1,2. However, an individual single-cell modality only reflects a snapshot of genetic features and partially depicts the peculiarity of cells, leading to characterization biases in complex biological systems2,3. Single-cell multi-omics (scMulti-omics) allows the quantification of multiple modalities simultaneously to fully capture the intricacy of complex molecular mechanisms and cellular heterogeneity. Such analyses advance various biological studies when paired with robust computational analysis methods4.

The existing tools for integrative analyses of scMulti-omics data, e.g., Seurat5, MOFA+6, Harmony7, and totalVI8, reliably predict cell types and states, remove batch effects, and reveal relationships or alignments among multiple modalities. However, most existing methods do not explicitly consider the topological information sharing among cells and modalities. Hence, they cannot effectively infer the active biological networks of diverse cell types simultaneously with cell clustering and have limited power to elucidate the response of these complex networks to external stimuli in specific cell types.

Recently, graph neural networks (GNN) have shown strength in learning low-dimensional representations of individual cells by propagating neighbor cell features and constructing cell-cell relations in a global cell graph9,10. For example, our in-house tool scGNN, a GNN model, has demonstrated superior cell clustering and gene imputation performance based on large-scale scRNA-seq data11. Furthermore, a heterogeneous graph with different types of nodes and edges has been widely used to model a multi-relational knowledge graph12. It provides a natural representation framework for integrating scMulti-omics data and learning the underlying cell-type-specific biological networks. Moreover, a recent development in the attention mechanism for modeling and integrating heterogeneous relationships has made deep learning models explainable and enabled the inference of cell-type-specific biological networks12,13.

In this work, we developed DeepMAPS (Deep learning-based Multi-omics Analysis Platform for Single-cell data), a heterogeneous graph transformer framework for cell-type-specific biological network inference from scMulti-omics data. This framework uses an advanced GNN model, i.e., heterogeneous graph transformer (HGT), which has the following advantages: (i) It formulates an all-in-one heterogeneous graph that includes cells and genes as nodes, and the relations among them as edges. (ii) The model captures both neighbor and global topological features among cells and genes to construct cell-cell relations and gene-gene relations simultaneously9,14,15,16. (iii) The attention mechanism in this HGT model enables the estimation of the importance of genes to specific cells, which can be used to discriminate gene contributions and enhances biological interpretability. (iv) This model is hypothesis-free and does not rely on the constraints of gene co-expressions, thus potentially inferring gene regulatory relations that other tools usually cannot find. It is noteworthy that DeepMAPS is implemented into a code-free, interactive, and non-programmatic interface, along with a Docker, to lessen the programming burden for scMulti-omics data.

Results

Overview of DeepMAPS

Overall, DeepMAPS is an end-to-end and hypotheses-free framework to infer cell-type-specific biological networks from scMulti-omics data. There are five major steps in the DeepMAPS framework (Fig. 1 and Methods). (i) Data are preprocessed by removing low-quality cells and lowly-expressed genes, and then different normalization methods are applied according to the specific data types. An integrated cell-gene matrix is generated to represent the combined activity of each gene in each cell. Different data integration methods are applied for different scMulti-omics data types5,6,7,8. (ii) A heterogeneous graph is built from the integrated matrix, including cells and genes as nodes and the existence of genes in cells as edges. (iii) An HGT model is built to jointly learn the low-dimensional embedding for cells and genes and generate an attention score to indicate the importance of a gene to a cell. (iv) Cell clustering and functional gene modules are predicted based on HGT-learned embeddings and attention scores. (v) Diverse biological networks, e.g., gene regulatory networks (GRN) and gene association networks, are inferred in each cell type.

a The overall framework of DeepMAPS. Five main steps were included in carrying out cell clustering and biological gene network inference from the input scMulti-omics data. b The graph autoencoder was inserted with a HGT model. The integrated cell-gene matrix was used to build a heterogeneous graph include all cells (green) and genes (purple) as nodes. The HGT model is trained on multiple subgraphs (50 subgraphs as an example) that cover nodes in the whole graph as many as possible. Each subgraph is used to train the model with 100 epochs; thus, the whole training process iterates 5,000 times. The trained model is then applied to the whole graph to learn and update the embeddings of each node. c An illustration of embedding update process of the target node in a single HGT layer. The red circle in the upper panel indicates the target node and the black circle indicates the source nodes. Arrows represents for the connection between a target node and source nodes. Colored rectangles represent for embeddings of different nodes. The zoom in detailed process in the bottom panel shows the massage passing process and attention mechanism. The final output of one HGT layer is an update of node embedding for all nodes. HGT heterogeneous graph transformer.

To learn joint representations of cells and genes, we first generate a cell-gene matrix integrating the information of the input scMulti-omics data. A heterogeneous graph with cell nodes and gene nodes is then constructed, wherein an unweighted cell-gene edge represents the existence of gene activity of a gene in a cell, and the initial embedding of each node is learned from the gene-cell integrated matrix via two-layer GNN graph autoencoders (Methods). Such a heterogeneous graph offers an opportunity to clearly represent and organically integrate scMulti-omics data so that biologically meaningful features can be learned synergistically. The entire heterogeneous graph is then sent to a graph autoencoder to learn the relations between the cells and genes and update the embedding of each node. Here, DeepMAPS adopts a heterogeneous multi-head attention mechanism to model the overall topological information (global relationships) and neighbor message passing (local relationships) on the heterogeneous graph. The heterogeneous graph representation learning provides a way to enable the embedding of cells and genes simultaneously using the transformer in DeepMAPS. The initial graph determines the path of message passing and how the attention scores can be calculated in DeepMAPS.

In each HGT layer, each node (either a cell or a gene) is considered a target, and its 1-hop neighbors as sources. DeepMAPS evaluates the importance of its neighbor nodes and the amount of information that can be passed to the target based on the synergy of node embedding (i.e., attention scores). As a result, cells and genes with highly positive correlated embeddings are more likely to pass messages within each other, thus maximizing the similarity and disagreement of the embeddings. To make the unsupervised training process feasible on a large heterogeneous graph, DeepMAPS is performed on 50 subgraphs sampled from the heterogeneous graph, covering a minimum of 30% of cells and genes to train for the shared parameters between different nodes, information which is later used for testing of the whole graph. As an important training outcome, an attention score is given to represent the importance of a gene to a cell. A high attention score for a gene to a cell implies that the gene is relatively important for defining cell identity and characterizing cell heterogeneity. This discrimination allows for the construction of reliable gene association networks in each cell cluster as the final output of DeepMAPS. We then build a Steiner Forest Problem (SFP) model17 to identify genes with higher attention scores and similar embedding features to a cell cluster. The gene-gene and gene-cell relations in the optimized solution of the SFP model mirror the embedding similarity of genes and the attention importance of genes to each cell cluster. A gene association network can be established by genes with the highest importance in characterizing the identity of that cell cluster based on their attention scores and embedding similarities, and these genes are considered to be cell-type-active.

DeepMAPS achieves superior performances in cell clustering and biological network inference from scMulti-omics data

We benchmarked the cell clustering performance of DeepMAPS on ten scMulti-omics datasets, including three multiple scRNA-seq datasets (R-bench-1, 2, and 3), three CITE-seq datasets (C-bench-1, 2, and 3), and four matched scRNA-seq and scATAC-seq (scRNA-ATAC-seq) datasets measured from the same cell (A-bench-1, 2, 3, and 4) (Supplementary Data 1). Specifically, the six R-bench and C-bench datasets have benchmark annotations provided in their original manuscripts, while the four A-bench datasets do not. These datasets cover the number of cells ranging from 3,009 to 32,029; an average read depth (considering scRNA-seq data only) ranging from 2,885 to 11,127; and a zero-expression rate (considering scRNA-seq data only) from 82 to 96% (Supplementary Data 1).

We compared DeepMAPS with four benchmarking tools (Seurat v3 and v45,18, MOFA + 6, TotalVI8, Harmony7, and GLUE19 (Methods)) in terms of the Average Silhouette Width (ASW), Calinski-Harabasz (CH), Davies-Bouldin Index (DBI), and Adjusted Rand Index (ARI) to evaluate cell clustering performance. For each dataset, we trained DeepMAPS on 36 parameter combinations, including the number of heads, learning rate, and the number of training epochs. To ensure fairness, each benchmarking tool was also tuned with different parameter combinations (Methods). DeepMAPS achieves the best performance comparing all benchmark tools in all test datasets in terms of ARI (for R-benches and C-benches) and ASW (for A-benches) (Fig. 2a, Supplementary Figs. 1–3, and Source Data 1-3). We also noticed that Seurat was the second-best performed tool, with small variances for different parameter selections in all benchmark datasets. We selected the default parameter per data type based on the performance of parameter combinations on the grid-search benchmarking. The parameter combination with the highest median ARI/ASW scores averaged in all benchmark datasets was considered as the default parameters for the corresponding data type.

a Benchmark cell clustering results of ten datasets in ARI for the three multiple scRNA-seq data and the three CITE-seq data with benchmark labels, and ASW for the four scRNA-ATAC-seq data without benchmark labels. Each box showcases the minimum, first quartile, median, third quartile, and maximum ARI or AWS results of a tool using different parameter settings (DeepMAPS: n = 96, Seurat: n = 16 for RNA-RNA and CITE-seq and 36 for RNA-ATAC, Harmony: n = 36, MOFA + : n = 36, TotalVI: n = 48, and GLUE: n = 72). Dots represent outliers. b Results comparison on five independent datasets. No repeated experiment was conducted. c Robustness test of DeepMAPS using the cell cluster leave-out method for the three independent test datasets with benchmarking cell labels. p-values were calculated based on two-tail t.test. Each box showcases the minimum, first quartile, median, third quartile, and maximum ARI results of a tool performed on different data subsets (R-test: n = 5, C-test: n = 20, and A-test-1: n = 5). Dots represent outliers. d–f UMAP comparison of R-test, C-test, and A-test-1 datasets between DeepMAPS and other tools using the original cell labels. Source data are provided as a Source Data file. ASW average Silhouette width, ARI adjusted rand index.

Additional benchmarking experiments were carried out to justify the selection of different integration methods in DeepMAPS. Specifically, for the analysis of scRNA-ATAC-seq data, we designed an integration method using gene velocity to balance the weight between gene expressions and chromatin accessibilities in characterizing cell activities and states (Methods). This integration process can ensure harmonizing datasets (especially for multiple scRNA-seq data) and generate an integrated matrix (with genes as rows and cells as columns) as the input for HGT. Our results showed that, for benchmark data 1 and 2 (A-bench-1 and −2), the velocity-based approach showed significantly (p-value <0.05) higher ASW scores than the weighted nearest neighbor (WNN) approach in Seurat v 4.0 on all grid-search parameter combinations (Supplementary Fig. 4 and Source Data 4). We reason that with the inclusion of velocity information, the modality weight between gene expression and chromatin accessibility that contribute to recognize cell types are better balanced (Supplementary Fig. 5). The comparison of modality weight of scATAC-seq in different cell clusters by using or without using the velocity-weighted balance method. In addition, we compared different clustering methods (i.e., Leiden, Louvain, and SLM) in DeepMAPS and compared the impact of clustering resolutions (i.e., 0.4, 0.8, 1.2, and 1.6) to cell clustering results. We found no significant differences among these clustering methods, and Louvain showed slightly better performance than the other two (Supplementary Fig. 6 and Source Data 5). Lastly, DeepMAPS achieved higher scores than other tools when selecting the same clustering resolution. We also found that, in most cases, higher resolution lower cell clustering prediction scores; therefore, we selected resolution at 0.4 as the default parameter in DeepMAPS (Supplementary Fig. 7 and Source Data 6–8).

We further independently tested our default parameter selection on five independent datasets, named R-test, C-test, A-test-1, −2, and −3, by comparing our results with the same benchmarking tools using their default parameters. For the three test datasets with benchmarking cell labels, DeepMAPS performs the best in terms of ARI score, while for the two scRNA-ATAC-seq datasets without cell labels, the benchmarking tools in the comparison achieve similar performance (Fig. 2b and Source Data 9). In order to evaluate the robustness of DeepMAPS, a leave-one-out test was performed on the three independent test datasets with benchmark labels (Fig. 2c and Source Data 10). We first removed all cells in a cell cluster based on benchmark labels and then applied DeepMAPS and other tools on the remaining cells. For each dataset, the leave-one-out results of DeepMAPS were better than the other tools with higher ARI scores, indicating that the message passing and attention mechanism used in DeepMAPS maintains cell-cell relations in a robust manner.

The cell clustering UMAP on the three independent datasets with benchmarking labels showcased that the latent representations obtained in DeepMAPS can better preserve the heterogeneity of scRNA-seq data (Fig. 2d–f). For the R-test dataset, all tools showed the ability to separate mesenchymal, leukocyte, and endothelial cells, but failed to separate urothelium basal cells and bladder cells. However, cells on the DeepMAPS UMAP are more compact, and the bladder cells (red dots) are grouped better than MOFA + and Seurat (Fig. 2d). For the C-test dataset, cells in the same cluster are more ordered and compact (e.g., the B cell cluster and NK cell cluster), while cells from different clusters are more apart from each other on DeepMAPS UMAP (e.g., CD8 cell clusters and CD4 cell clusters). (Fig. 2e). For the A-test-1 dataset, DeepMAPS was the only tool that accurately separated each cell type. In contrast, Seurat and MOFA + mistakenly divided the PDX1 or PDX2 population into two clusters and included more mismatches (Fig. 2f).

DeepMAPS can infer statistically significant and biologically meaningful gene association networks from scMulti-omics data

We evaluated the two kinds of biological networks that DeepMAPS can infer, gene association network and GRN, in terms of centrality scores and functional enrichment. For the R-test dataset (Fig. 3a) and C-test dataset (Fig. 3b), we used two centrality scores, closeness centrality (CC) and eigenvector centrality (EC), that have been used in previous single-cell gene association network evaluations20, to compare the identified gene association networks from all the tools in this comparison. CC reflects the average connectivity of a node to all others in a network, and EC reflects the importance of a node based on its connected nodes. Both CC and EC can interpret node’s influence in identifying genes that may play more critical roles in the network. A gene association network with higher node centrality indicates that the detected genes are more likely to be involved in critical and functional biological systems. We also constructed gene co-expression networks as a background using all genes in a dataset by calculating Pearson’s correlation coefficients of gene expression in a cell cluster. p-value = 0.05 was set as the edge cutoff. We compared cell-type-active gene association networks generated in DeepMAPS with those generated in IRIS315 and the background co-expression networks. The average CC and EC of networks constructed by DeepMAPS in R-test and C-test datasets showed significantly higher scores than IRIS3 and the background co-expression networks (Source Data 11). We reason that the gene association network generated in DeepMAPS is not only co-expressed but also of great attention impact on cells; thus, genes in the network tend to be more important to the cell type.

a, b Closeness centrality (CC) and eigenvector centrality (EC) were used to indicate the compactness and importance of genes to the network. We compared our results with IRIS3 and a background network using all genes for the R-test dataset (n = 5) a and C-test dataset (n = 14) b. p-values were calculated using a two-tail t-test. c Comparison of the number of unique TFs in GRNs that showed significantly enriched biological functions in three public databases. Each box contains the results of six scRNA-ATAC-seq datasets (n = 6). d Comparison of the number of cell-type-specific regulons in GRNs significantly enriched in only one biological function/pathway in the three public databases (n = 6). e The F1 score comparisons of regulons enriched to only one function/pathway using three databases (n = 6). The mean value of precision and recall scores of the selected six scRNA-ATAC-seq datasets were max-min scaled and shown in the heatmap with darker blue indicating high values and lighter blue indicating low values. Source data are provided as a Source Data file. Each box in Fig. 3 showcases the minimum, first quartile, median, third quartile, and maximum score of the corresponding criteria. CC closeness centrality, EC eigenvector centrality, CTSR cell-type-specific regulon.

To evaluate whether DeepMAPS can identify biologically meaningful GRNs in specific cell types, we performed enrichment tests on basic gene regulatory modules (i.e., regulons14), with three public functional databases, Reactome21, DoRothEA22, and TRRUST v223. To avoid any bias in comparison, we compared cell-type-specific GRNs inferred from DeepMAPS with (i) IRIS3 and SCENIC14 on the scRNA-seq matrix, (ii) IRIS3 and SCENIC on a gene-cell matrix recording the gene activity scores (GAS) calculated in DeepMAPS based on the velocity-based integration method, (iii) MAESTRO24 on scATAC-seq matrix, and (iv) MAESTRO on original scRNA-seq and scATAC-seq matrix. The six datasets collected from human tissue were used (i.e., A-test-1, A-bench-2, A-bench-3, A-bench-4, A-test-1, A-test-2). We first showcased that the GRNs identified in DeepMAPS included more unique transcription factor (TF) regulations than the other tools, except for the enrichment to the DoRothEA database (Fig. 3c and Source Data 12). We considered that a highly cell-type-specific regulon (CTSR) might represent only one significant enriched functionality; alternatively, a generic regulon might improperly contain genes involved in several pathways. Therefore, we compared the number of CTSRs enriched to one function/pathway across different tools. DeepMAPS outperformed (p-value<0.05) other tools on most of the six scRNA-ATAC-seq datasets in terms of the number of regulons that enrich only one function/pathway and the enrichment F1 scores (Fig. 3d, e and Source Data 12). For the F1 score of the enrichment test to the TRRUST v2 database, DeepMAPS (median F1 score is 0.026) was slightly lower than IRIS3 using the GAS matrix (median F1 score is 0.031). We also noticed that all tools did not achieve good enrichments in the TRRUST v2 database mainly due to the small number of genes (on average, 10 genes are regulated by one TF; 795 TFs in total). SCENIC also showed competitive scaled precision scores (scaled mean: 0.47 for Reactome, 0.66 for DoRothEA, and 0.61 for TRRUST v2), while achieving lower scaled recall scores, making the F1 scores smaller than DeepMAPS for most datasets. IRIS3 and SCENIC performed on the GAS matrix showed better enrichment results than using scRNA-seq data only, indicating that integrating information from scRNA-ATAC-seq data is more useful for GRN inference than using scRNA-seq data alone.

DeepMAPS accurately identifies cell types and infers cell-cell communication in PBMC and lung tumor immune CITE-seq data

We present a case study that applies DeepMAPS to a published mixed peripheral blood mononuclear cells (PBMC) and lung tumor leukocytes CITE-seq dataset (10× Genomics online resource, Supplementary Data 1) to demonstrate capacity in modeling scMulti-omics in characterizing cell identities. The dataset includes RNAs and proteins measured on 3485 cells. DeepMAPS identified 13 cell clusters, including four CD4+ T cell groups (naive, central memory (CM), tissue-resident memory (TRM), and regulatory (Treg)), two CD8+ T cell groups (CM and TRM), a natural killer cell group, a memory B cell group, a plasma cell group, two monocyte groups, one tumor-associated macrophage (TAM) group, and a dendritic cell (DC) group. We annotated each cluster by visualizing the expression levels of curated maker genes and proteins (Fig. 4a, b and Supplementary Data 2). Compared to cell types identified using only proteins or RNA, we isolated or accurately annotated cell populations that could not be characterized using the individual modality analysis. For example, the DC cluster was only successfully identified using the integrated protein and RNA. By combining signals captured from both RNA and proteins, DeepMAPS successfully identified biologically reasonable and meaningful cell types in the CITE-seq data.

a UMAPs for DeepMAPS cell clustering results from integrated RNA and protein data, protein data only, and RNA data only. Cell clusters were annotated based on curated marker proteins and genes. b Heatmap of curated marker proteins and genes that determine the cell clustering and annotation. c Heatmap of the Spearman correlation comparison of top differentially expressed genes and proteins in plasma cells and memory B cells. d UMAP is colored by the 51st embedding, indicating distinct embedding representations in plasma cells and memory B cells. e Expression of top differentially expressed genes and proteins in c as a function of the 51st embedding to observe the pattern relations between plasma cells and memory B cells. Each line represents a gene/protein, colored by cell types. For each gene, a line was drawn using a loess smoothing function based on the corresponding embedding and scaled gene expression in a cell. f–h Similar visualization was conducted for the 56th embedding to compare EM CD8+ T cells and TRM CD8+ T cells c–e. i Two signaling pathways, NECTIN and ALCAM, are shown to indicate the predicted cell–cell communications between two cell clusters. A link between a filled circle (resource cluster with highly expressed ligand coding genes) and an unfilled circle (target cluster with highly expressed receptor coding genes) indicates the potential cell-cell communication of a signaling pathway. Circle colors represent different cell clusters, and the size represents the number of cells. The two monocyte groups were merged. TRM tissue-resident memory, CM central memory, TAM tumor-associated macrophage, HGT heterogeneous graph transformer.

We then compared the modality correlation between the two cell types. We used the top differentially expressed genes and proteins between memory B cells and plasma cells, and performed hierarchical clustering of the correlation matrix. The result clearly stratified these features into two anticorrelated modules: one associated with memory B cells and the other with plasma cells (Fig. 4c). Furthermore, we found that the features in the two modules significantly correlated with the axis of maturation captured by our HGT embeddings (Supplementary Fig. 8 and Supplementary Data 3). For example, one HGT embedding (the 51st) showed distinctive differences between plasma cells and memory B cells (Fig. 4d, e). Similar findings were also observed when comparing EM CD8+ T cells with TRM CD8+ T cells (Fig. 4f). Nevertheless, it is possible to identify a representative HGT embedding (56th) that maintains embedding signals for a defined separation of the two groups (Fig. 4g, h). These results point to any two cell clusters consisting of coordinated activation and repression of multiple genes and proteins, leading to a gradual transition in cell state that can be captured by a specific dimension of the DeepMAPS latent HGT space. On the other hand, we generated the gene-associated networks with genes showing high attention scores for EM CD8+ T cells, TRM CD8+ T cells, memory B cells, and plasma cells and observed diverse patterns (Supplementary Fig. 9).

Based on the cell types and raw data of gene and protein expressions, we inferred cell–cell communications and constructed communication networks among different cell types within multiple signaling pathways using CellChat25 (Fig. 4i). For example, we observed a CD6-ALCAM signaling pathway existing between DC (source) and TRM CD4+ T cells (target) in the lung cancer tumor microenvironment (TME). Previous studies have shown that ALCAM on antigen-presenting DCs interacts with CD6 on the T cell surface and contributes to T cell activation and proliferation26,27,28. As another example, we identified the involvement of the NECTIN-TIGIT signaling pathway during the interaction between the TAM (source) and TRM CD8+ T cells (target), which is supported by a previous report that NECTIN (CD155) expressed on TAM could be immunosuppressive when interacting with surface receptors, TIGIT, on CD8+ T cells in the lung cancer TME29,30.

DeepMAPS identifies specific GRNs in diffuse small lymphocytic lymphoma scRNA-seq and scATAC-seq data

To further extend the power of DeepMAPS to GRN inference, we used a single-cell Multiome ATAC + Gene expression dataset available on the 10× Genomics website (10× Genomics online resource). The raw data is derived from 14,566 cells of flash-frozen intra-abdominal lymph node tumor from a patient diagnosed with diffuse small lymphocytic lymphoma (DSLL) of the lymph node lymph. We integrate gene expression and chromatin accessibility by balancing the weight of each modality of a gene in a cell based on RNA velocity (Fig. 5a and Method). To build TF-gene linkages, we considered gene expression, gene accessibility, TF-motif binding affinity, peak-to-gene distance, and TF-coding gene expression. Genes found to be regulated by the same TF in a cell cluster are grouped as a regulon. We considered regulons with higher centrality scores to have more significant influences on the characterization of the cell cluster. Regulons regulated by the same TF across different cell clusters are compared for differential regulon activities. Those with significantly higher regulon activity scores (RAS) are considered as the cell-type-specific regulons in the cell cluster.

a Conceptual illustration of DeepMAPS analysis of scRNA-ATAC-seq data. Modalities are first integrated based on a velocity-weighted balance. The integrated GAS matrix was then used to build a heterogeneous graph as input into the HGT framework. The cell cluster and gene modules with high attention scores were then used for building TF-gene linkages and determining regulons in each cell cluster. b The UMAP shows the clustering results of DeepMAPS. Cell clusters were manually annotated based on curated marker genes. c The observed and extrapolated future states (arrows) based on the RNA velocity of the normal B cell and the two DSLL states are shown (top panel). Velocity-based trajectory analysis shows the pseudotime from the top to the bottom right (bottom panel). d Selected 20 TF in each of the three clusters, representing the top 20 regulons with the highest centrality scores. Colors represent regulons uniquely identified in each cluster or shared between different clusters. e Regulons in DSLL state-1 showed a significant difference in regulon activity compared to the other clusters. Motif shape and number of regulated genes are also shown. f Violin plots of regulon activities of the four regulons compared between the three clusters. g The downstream-regulated genes of JUN (the most differentially active regulon in DSLL state-1) in the three clusters. h An illustration of the BAFF signaling pathway identified from GAS-based cell-cell communication prediction using CellChat. The BAFF signaling pathway was found to exist between macrophage and both DSLL states. It further activates the JUN regulon and enables the transcription of genes like CDK6. Figure created with BioRender.com. i The ATAC peak, RNA expression, and GAS level of TNFRSF13B (the coding gene of TACI, the receptor in the BAFF signaling pathway). Source data are provided as a Source Data file.

DeepMAPS identified 11 cell clusters in the DSLL data. All clusters were manually annotated based on curated gene markers (Fig. 5b and Supplementary Data 4). Two DSLL-like cell clusters (DSLL state-1 and state-2) were observed. The RNA velocity-based pseudotime analysis performed on the three B cell clusters (normal B cell and two DSLL states) assumed that the two DSLL states were derived from normal B cells, and state-1 was derived earlier than state-2, although the two states seemed to be partially mixed (Fig. 5c). We further selected the top 20 TFs with the highest regulon centrality scores in each of the three cell clusters (Fig. 5d and Source Data 13). Interestingly, these TFs showed distinctions between the normal and the two DSLL states and inferred variant regulatory patterns within the two DSLL states. For regulons shared by all three B cell clusters, EGR1, MEF2B, and FOS were transcriptionally active in both normal B and DSLL cells and responsible for regulating B cell development, proliferation, and germinal center formation31,32,33,34. E2F6, ELF3, and KLF16 were identified as shared only in the two DSLL states, with reported roles in tumorigenesis35,36,37,38,39,40. Further, JUN, MAFK, and MAFG, which encode the compartments of the activating protein-1 (AP-1),34,41,42 were found to be active in DSLL state-1 while NFKB1, coding for a subunit of the NF-κB protein complex43,44, was found to be active in DSLL state-2.

We constructed a GRN consisting of the four cell-type-specific regulons (JUN, KLF16, GATA1, and FOS) (Fig. 5e and Supplementary Fig. 10) in DSLL state-1 with RAS that is significantly higher than normal B cells and DSLL state-2 (Fig. 5f). KLF16 reportedly promotes the proliferation of both prostate39 and gastric cancer cells40. FOS and JUN are transcription factors in the AP-1 family, regulating the oncogenesis of multiple types of lymphomas34,41,42,45, and GATA1 is essential for hematopoiesis, the dysregulation implicated in multiple hematologic disorders, and malignancies46,47. Distinct regulatory patterns were also observed when we zoomed in on a single regulon (Fig. 5g and Supplementary Figs. 11-12). As the most active regulon in DSLL state-1, JUN was found to regulate five unique downstream genes and 12 genes shared with DSLL state-2. Downstream genes, including CDK633,34, IGF2R48, and RUNX149, are critical for cell proliferation, survival, and development functions in DSLL.

Moreover, we further built connections between upstream cell-cell communication signaling pathways and downstream regulatory mechanisms in DSLL cells. We identified a cell-cell communication between macrophage and the two DSLL states via the B cell activation factor (BAFF) signaling pathway, based on the integrated GAS matrix using CellChat25, which includes BAFF as the ligand on macrophage cells and TACI (transmembrane activator and calcium-modulator and cyclophilin ligand interactor) as the receptor on DSLL cells (Fig. 5h). BAFF signaling is critical to the survival and maturation of normal B cells50,51, while aberrations contribute to the resistance of malignant B cells to apoptosis52,53. We observed that the expression of the TACI coding gene, TNFRSF13B, was explicitly higher in the two DSLL states, while the corresponding chromatin accessibility maintained high peaks in state-1 (Fig. 5i). Upon engagement with its ligand, TACI has been reported to transduce the signal and eventually activate the AP-154,55 and NF-κB56,57 transcriptional complexes for downstream signaling in B cells. JUN (a subunit of AP-1) was identified as the most specific and key regulator in state-1 responsible for cell proliferation and regulating downstream oncogenes, such as CDK6, that has been reported to promote the proliferation of cancer cells in multiple types of DSLLs as well as other hematological malignancies58,59,60. It is clear that BAFF signaling first appears in DSLL state-1 and triggers the activation of the JUN regulatory mechanism, leading to a high regulon activity of JUN. The JUN regulon accelerates the proliferation and oncogenesis explicitly in DSLL, leading to a more terminal differential stage of DSLL (state-2). As a result, state-1 includes cells undergoing rapid cell proliferation and differentiation, transitioning from normal B cells to matured DSLL. In short, DeepMAPS can construct GRNs and identify cell-type-specific regulatory patterns to offer a better understanding of cell states and developmental orders in diseased subpopulations.

DeepMAPS provides a multi-functional and user-friendly web portal for analyzing scMulti-omics data

Due to the complexity of single-cell sequencing data, more webservers and dockers have been developed in the past three years61,62,63,64,65,66,67,68,69,70,71,72,73 (Supplementary Data 5). However, most of these tools only provide minimal functions such as cell clustering and differential gene analysis. They do not support the joint analysis of scMulti-omics data and especially lack sufficient support for biological network inference. On the other hand, we recorded the running time of DeepMAPS and benchmark tools on different datasets with cell numbers ranging from 1000 to 160,000 (Supplementary Data 6). The deep learning models (DeepMAPS and TotalVI) have longer running time than Seurat and MOFA + . To these ends, we provided a code-free, interactive, and non-programmatic interface to lessen the programming burden for scMulti-omics data (Fig. 6a). The webserver supports the analysis of multiple RNA-seq data, CITE-seq data, and scRNA-ATAC-seq data using DeepMAPS (Fig. 6b). Some other methods, e.g., Seurat, are also incorporated as an alternative approach for the users’ convenience. Three major steps—data preprocessing, cell clustering and annotation, and network construction—are included in the server. In addition, the DeepMAPS server supports real-time computing and interactive graph representations. Users may register for an account to have their own workspace to store and share analytical results. Other than the advances mentioned, the DeepMAPS webserver highlights an additional function for the elucidation of complex networks in response to external stimuli in specific cell types. The user can upload a metadata file with phenotype information (e.g., cells with treatment and without treatment), select, and re-label the corresponding cells (e.g., CD8+ T cells with treatment and CD8+ T cells without treatment). In this way, DeepMAPS will predict the treatment-related networks in CD8+ T cells. Examples are given in the online tutorial at https://bmblx.bmi.osumc.edu/tutorial.

Discussion

DeepMAPS is a deep-learning framework that implements heterogeneous graph representation learning and a graph transformer in studying biological networks from scMulti-omics data. By building a heterogeneous graph containing both cells and genes, DeepMAPS identifies their joint embedding simultaneously and enables the inference of cell-type-specific biological networks along with cell types in an intact framework. Furthermore, the application of a heterogeneous graph transformer models the cell-gene relation in an interpretable uniform multi-relation. In such a way, the training and learning process in a graph can be largely shortened to consider cell impacts from a further distance.

By jointly analyzing gene expression and protein abundance, DeepMAPS accurately identified and annotated 13 cell types in a mixed CITE-seq data of PBMC and lung tumor leukocytes based on curated markers that cannot be fully elucidated using a single modality. We have also proved that the embedding features identified in DeepMAPS capture statistically significant signals and amplify them when the original signals are noisy. Additionally, we identified biologically meaningful cell-cell communication pathways between DC and TRM CD4+ T cells based on the gene association network inferred in the two clusters. For scRNA-ATAC-seq, we employed an RNA velocity-based method to dynamically integrate gene expressions and chromatin accessibility that enhanced the prediction of cell clusters. Using this method, we identified distinct gene regulatory patterns among normal B cells and two DSLL development states. We further elucidated the deep biological connections between cell-cell communications and the downstream GRNs, which helped characterize and define DSLL states. The identified TFs and genes can be potential markers for further validation and immuno-therapeutical targets in DSLL treatment.

While there are advantages and improved performances for analyzing scMulti-omics data, there is still room to improve the power of DeepMAPS further. First, the computational efficiency for super-large datasets (e.g., more than 1 million cells) might be a practical issue considering the complexity of the heterogeneous graph representation (which may contain billions of edges). Moreover, DeepMAPS is recommended to be run on GPUs, which leads to a potential problem of reproducibility. Different GPU models have different floating-point numbers that may influence the precision of loss functions during the training process. For different GPU models, DeepMAPS may generate slightly different cell clustering and network results. Lastly, the current version of DeepMAPS is based on a bipartite heterogeneous graph with genes and cells. Separate preprocessing and integration steps are required to transfer different modalities to genes for integration into a cell-gene matrix. To fully achieve an end-to-end framework for scMulti-omics analysis, the bipartite graph can be extended to a multipartite graph, where different modalities can be included as disjoint node types (e.g., genes, proteins, or peak regions). Such a multipartite heterogeneous graph can also include knowledge-based biological information, such as known molecular regulations and more than two modalities in one graph. However, by including more node types, the computational burden will be increased geometrically, which requires a dedicated discovery of model and parameter optimization in the future.

In summary, we evaluated DeepMAPS as a pioneer study for the integrative analysis of scMulti-omics data and cell-type-specific biological network inference. It will likely provide different visions of deep learning deployment in single-cell biology. With the development and maintenance of the DeepMAPS webserver, our long-term goal is to create a deep learning-based eco-community for archiving, analyzing, visualizing, and disseminating AI-ready scMulti-omics data.

Abstract