- Method
- Open access
- Published:

# MNMST: topology of cell networks leverages identification of spatial domains from spatial transcriptomics data

*Genome Biology*
**volume 25**, Article number: 133 (2024)

## Abstract

Advances in spatial transcriptomics provide an unprecedented opportunity to reveal the structure and function of biology systems. However, current algorithms fail to address the heterogeneity and interpretability of spatial transcriptomics data. Here, we present a multi-layer network model for identifying spatial domains in spatial transcriptomics data with joint learning. We demonstrate that spatial domains can be precisely characterized and discriminated by the topological structure of cell networks, facilitating identification and interpretability of spatial domains, which outperforms state-of-the-art baselines. Furthermore, we prove that network model offers an effective and efficient strategy for integrative analysis of spatial transcriptomics data from various platforms.

## Background

Complex tissues are typically composed of many cell types (sub-populations) located at various positions to execute the corresponding biological functions [1,2,3]. Thus, it is of great significance to simultaneously exploit the expression and spatial information of cells, which is the foundation for understanding the underlying mechanisms of biology systems [4]. For example, the combination of expression and spatial information of cells provides an effective strategy to precisely characterize the heterogeneity of tumor micro-environment, which is critical for diagnosis and therapy of cancers [5]. Fortunately, advances in biological techniques ensure the generation of expression profiles of cells while retaining spatial relative context in situ, known as spatial transcriptomics, which presents unprecedented opportunities to explore the organization and function of tissues in a spatial context. On the basis of principles of preserving spatial information, current technologies are broadly classified into two categories, i.e., imaging- and next-generation sequencing (NGS)-based methods [6]. The typical imaging-based approaches include fluorescence in situ hybridization (FISH) [7] and its variants [8,9,10], which are criticized for their limited capacity to detect RNA transcripts. To address this issue, NGS-based methods employ spatial barcoding and NGS, including Slide-seq [11], Stereo-seq [12], and 10 \(\times\) Visium [13]. The accumulated spatial transcriptomics data provide a great opportunity to investigate the functions and cellular structure of tissues by exploiting interesting patterns and features that cannot be identified from other data.

Spatial domains, that are regions where cells inside are coherent with similar expression profiles and spatial proximity, are one of the typical patterns in spatial transcriptomics data. They are the prerequisites for downstream analysis, such as tracking disease progression [14] and tissue development [8]. According to strategies of dissection, available methods for the identification of spatial domains are divided into two categories, i.e., biological experiments and computational approaches. Specifically, former methods employ anatomy to obtain the spatial information of cells and manually select spatial domains. This approach is precise and reliable given that expertise knowledge is seamlessly incorporated. However, these methods are criticized for efficiency because manual selection mode poses a great challenge in terms of time and finance. Computational methods, which identify spatial domains with machine learning, provide an alternative to overcome limitations of biological experiment-based methods.

Specifically, unsupervised clustering is a widely adopted strategy for identifying spatial domains in spatial transcriptomics data. Current computational methods for spatial domain identification are roughly categorized into two classes, i.e., non-spatial and spatial clustering algorithms. Non-spatial clustering methods, such as DRjCC [15], SCANPY [16], K-means, and Louvain [17], are deliberately designed for transcriptomics data without spatial information. These methods are criticized for their undesirable performance because identified spatial domains deviate from tissue sections, highlighting the importance of spatial information for the analysis of spatial transcriptomics data.

However, spatial clustering integrates gene expression and spatial location to address the aforementioned issues. This immediate purpose is to account for spatial dependency of gene expression, which results in better matching spatial location. In comparison to non-spatial clustering, spatial clustering significantly enhances the performance of algorithms, which indicates that spatial information plays an indispensable role in the identification of spatial domains. The greatest difference among these algorithms lies in how to obtain and integrate spatial and transcriptional features of cells. For example, BayesSpace [18] and SpatialPCA [19] employ Hidden-Markov random field to enforce physically proximal cells belonging to the same domains by assigning higher probability to proximal cells. However, these algorithms achieve an excellent performance if and only if cell types are well separated, that is the division boundary is very clear. Actually, the boundaries of spatial domains are unclear because of the hierarchical structure of cell types.

To exploit the latent structure of cells, SpaGCN [20] and STAGATE [21] adopt graph neural networks to learn the topological structure of cells, which can be concatenated with transcriptional features to facilitate identification of spatial domains. These algorithms differ greatly in strategies used to learn graph features of cells for characterizing spatial information. For example, SEDR [22] employs auto-encoder to learn spatial embedding, whereas CCST [23] addresses complex global cell interactions across tissues. BANKSY [24] simultaneously considers transcriptome of cells and their local neighbors, and DeepST [25] integrates transcriptomics and morphological features of cells with data augmentation, which provides a compatible strategy for spatial domain. SPARK-X [26] and SpatialDE [27] detect spatial variable genes to link spatial domains with biological functions for further integrating knowledge from multiple resources.

Although great efforts have been devoted to identifying spatial domains, many urgent but unaddressed problems exist. Current algorithms directly integrate the spatial and expression data. However, they are criticized for failing to remove heterogeneity of spatial transcriptomics data, which results in an undesirable performance. Second, the transcriptional and spatial features of cells are independently learned, where the relations among spatial and expression information are neglected, which leads to a failure to precisely model spatial domains. Third, available algorithms are devoted to characterize and identify spatial domains with deep features of cells, reducing interpretability of cell features. Finally, vast majority of algorithms are designed for specific platforms, which hinders the application and transplanting to others. Thus, precisely modeling and identifying spatial domains from spatial transcriptomics data is still challenging.

We propose a multi-layer network model for spatial transcriptomics data (MNMST) to accurately identify spatial domains with joint learning for addressing the aforementioned issues. As shown in Fig. 1, MNMST consists of three major components: multi-layer network construction, spatial and expression feature learning, and spatial domain identification. Specifically, MNMST first constructs a network for spatial and expression data to overcome the heterogeneity of spatial transcriptomics data. This procedure results in homogeneous multi-layer networks. Subsequently, it transforms spatial domain identification into multi-layer network clustering problem. Then, it employs joint learning to integrate constructed multi-layer networks by projecting cells into a subspace, where the low-dimensional and compatible features of cells are learned. To enhance the interpretability of features, MNMST employs self-representation learning (SRL) to construct an affinity graph of cells with cell features. Experimental results demonstrate that the proposed multi-layer network model not only outperforms state-of-the-art algorithms on identification of spatial domains, but also covers all spatial transcriptomics platforms, including 10 \(\times\) Visium, osmFISH, STARMap, Stereo-seq, and Slide-seq V2. These results further prove that MNMST achieves the best performance for spatial domains identification, demonstrating that MNMST is promising for analyzing spatial transcriptomics data and understanding of tissue architecture and functions.

## Results

### Overview of the multi-layer network model for spatial domain identification

To smooth the reading and understanding of this study, we first describe the rationale of MNMST in this section (more technical description can be referred in the “Methods” section). For the sake of convenience, we use cells to represent measurement units in spatial transcriptomics, which can be freely replaced to spots in various platforms, such as 10\(\times\) Visium.

For any spatial transcriptomics data, the traditional clustering approaches [16, 28, 29] first perform feature learning for spatial and expression information of cells and then identify spatial domains by clustering the low-dimensional features of cells. However, these algorithms are criticized for the undesirable performance because features of cells fail to model and characterize the intrinsic and indirect relations among cells. To overcome this limitation, network-based algorithms, such as DeepST [25], SpaGCN [20], and STAGATE [21], first construct cell graph(s) by integrating spatial, expression information, etc. And, these algorithms then learn features by exploiting topological structure of the constructed cell networks, thereby dramatically enhancing performance of spatial domain identification.

However, these network-based algorithms still have several typical limitations that urgently need to be addressed. First of all, current algorithms construct cell network by calculating distances of cell pairs with concatenated cells features of spatial, expression, and morphological information, which are criticized for ignoring heterogeneity of spatial transcriptomics data, reducing quality and reliability of cell networks. Second, even though many alternatives are available for measuring distances among cells, it is difficult to select an appropriate manner to precisely characterize distances of cells (how measurement of distance effects performance of algorithms is studied in the “Systematic investigation of parameters of multi-layer network model” section). Finally, the learned low-dimensional features of cells only preserve the topological structure of cell networks, which neglect relations among spatial and expression information, failing to fully characterize the structure of spatial domains.

To address these issues, we propose a multi-layer network model (MNMST) to characterize and identify spatial domains in spatial transcriptomics data by integrating gene expression and spatial information of cells (Fig. 1). Specifically, to remove heterogeneity of spatial transcriptomics data, MNMST constructs two cell networks for spatial and expression data, respectively. In this case, MNMST converts the heterogeneous data integration problem into the homogeneous graphs clustering problem. Furthermore, MNMST learns features of cells by jointly exploiting topological structure of cell spatial and expression network, where the relation between spatial and expression information is implicitly learned, thereby improving quality of features of cells. In details, as shown in Fig. 1, MNMST consists of three major components, i.e., cell multi-layer network construction, cell feature learning, and spatial domain identification. On multi-layer network construction issue, MNMST first employs Euclidean distance to construct cell spatial network by exploiting indirect topological structure. Then, it automatically learns cell expression networks by utilizing self-representation learning (SRL) with local preservation constraint, where cell spatial information is incorporated to augment expression profiles of cells (Fig. 1 (A); see the “Methods” section). In this case, MNMST not only avoids selecting measurements of distances for cells, but also enhances quality of constructed graph.

On the cell feature learning issue, current algorithms independently learn spatial and expression features of cells, which ignores the latent relations among them. MNMST simultaneously learns spatial and expression features of cells to overcome this limitation. It does so by jointly decomposing cell multi-layer networks with nonnegative matrix factorization, which ensures the compatibility and quality of features (Fig. 1 (B)). MNMST automatically learns an affinity graph of cells with low-rank and sparse constraint to perform spatial domain identification. This way improves interpretability of spatial domains (Fig. 1 (B)). On the basis of affinity graph, MNMST identifies spatial domains and downstream analysis (Fig. 1 (C)).

Notice that, as far as we know, MNMST is the first multi-layer network model for spatial domain identification, which provides an alternative for modeling and analyzing spatial transcriptomics data.

### Benchmarking MNMST against feature-based state-of-the-art baselines

Human dorsolateral prefrontal cortex (DLPFC) data [30] consist of manually annotated cortical layers (L1–L6) and white matter (WM) of 12 slices with gene markers and cytoarchitecture. These data are generated with 10 \(\times\) Visium platform, serving as a publicly available benchmarking dataset for spatial domain identification (Fig. 2A). Seven feature-based state-of-the-art algorithms, namely SCANPY (with Leiden clustering) [16], Giotto [29], stLearn [28], BayesSpace [18], SpaGCN [20], DeepST [25], and SEDR [22], are selected as baselines to fully validate performance of MNMST. Notably, SCANPY is non-spatial, and others are spatial-based methods. Adjusted rand index (ARI) [31] is selected to measure the performance of algorithms (details in the “Methods” section).

Figure 2A visualizes slice 151675 of DLPFC with seven layers labeled with different colors, and MNMST is much more precise to identify these layers than baselines (Fig. 2B). Specifically, ARI of MNMST is 0.631, whereas that is 0.224 for Giotto, 0.263 for stLearn, 0.524 for BayesSpace, 0.461 for SpaGCN, and 0.521 for DeepST, respectively. Interestingly, SpaGCN, DeepST, and MNMST dramatically outperform others, demonstrating that topological structure of cell networks is also critical for characterizing spatial domains. However, DeepST and SpaGCN independently learn spatial and expression features of cells and explore relations of cell features with post-processing techniques. However, they ignore the latent relations among features. MNSMT jointly learns features of cells by decomposing matrices associated with multi-layer networks to overcome this limitation, where the latent relations of features are implicitly explored, thereby enhancing quality of features. Notice that MNMST is the only algorithm to discriminate L1 and L2 cortical layers, which cannot be delineated by others (Fig. 2B).

Furthermore, the proposed algorithm is also superior to baselines on spatial domain identification for other slices of DLPFC (Additional file 1: Figs. S1-S2), and Fig. 2C summarizes performance of various algorithms for DLPFC data, where *x*-axis denotes ARI to measure similarity of the predicted spatial layers and manually annotated layers. These results demonstrate that spatial algorithms significantly outperform non-spatial methods, showing that spatial information is critical need for spatial domain identification. In details, ARI of SCANPY is 0.208 ± 0.049 (median ± standard deviation), whereas ARI of the worst spatial algorithm is 0.250 ± 0.078. Specifically, ARI of MNMST is 0.553 ± 0.073, whereas ARI of the best baseline DeepST is 0.515 ± 0.001. Moreover, MNMST obtains the best performance in slice 151670 (Additional file 1: Fig. S1, ARI = 0.750), 151674 (Additional file 1: Fig. S2, ARI = 0.610), 151675 (Additional file 1: Fig. S2, ARI = 0.631), and 151676 (Additional file 1: Fig. S2, ARI = 0.588). These results prove that MNMST is superior to DeepST and SpaGCN, showing that the proposed multi-layer network model is much precise than current network ones. To investigate why multi-layer network model is more discriminative than baselines for spatial domain identification, PAGA [32] is employed to infer relations of identified spatial domains, where organization of various cortical layers is derived. Figure 2D illustrates trajectory of layers, where spatial domains identified by MNMST are well discriminated, and those by SCANPY are mixed, demonstrating that network model precisely exploits intrinsic structure of spatial domains.

Strikingly, we find that spatial domains can be clearly characterized by topology of cell affinity graph learned by MNMST. Figure 2E (left) visualizes topology structure of sub-graphs included by L1 and L2 layer, where connectivity is strong within layers, and weak across layers, corresponding to clusters of graphs. Four topological indexes, namely edge weight, degree (sum of weights on adjacent edges to a vertex), eigenvector centrality (importance of vertices with the largest eigenvector), and closeness (number of edges through a vertex), are selected to validate the difference between L1 and L2 layers for investigating the relation of spatial domains and topological structure of the affinity graph (right panel of Fig. 2E). Specifically, weights on edges across layers are much lower than those inside layers (*p *= 4.1E−10, Mann-Whitney *U* test). Furthermore, degrees, centrality, and closeness of cells in L1 and L2 layers significantly differ (degree: *p *= 7.3E−18, centrality: *p *= 3.9E−4, closeness: *p *= 4.0E−2, Mann-Whitney *U* test). These results demonstrate that the topological structure of cell networks is highly associated with spatial domains, providing an alternative for analyzing spatial transcriptomics data.

Next, spatial transcriptomics data of mouse brain tissue generated by 10 \(\times\) Visium platform is also adopted to evaluate the performance of MNMST. Brain anatomical references from the Allen Mouse Brain Atlas [33] are selected as baselines. MNMST precisely detects the cornu ammonis in mouse brain, and cerebellar cortex region in the sagittal posterior (Fig. 2F), which matches reference annotations [33]. Spatial domains identified by various algorithms for mouse brain tissue, and Silhouette Coefficient (SC) and Davies-Bouldin (DB) scores of them demonstrate that stLearn, DeepST, and MNMST are comparable (Additional file 1: Fig. S3). stLearn integrates spatial, transcriptional, and morphological information to comprehensively characterize structure of spatial domains in brain, whereas MNMST and DeepST fully exploit topological structure of cell networks for spatial domains. However, boundary of spatial domains identified by stLearn is non-smooth because it fails to fully exploit relations among spatial, transcriptional, and morphological information. These results show that the proposed multi-layer network model is also promising for characterizing complex spatial structure in mouse brain slices.

### Systematic investigation of parameters of multi-layer network model

MNMST learns features of cells with matrix factorization for multi-layer networks, implying that the theoretical complexity of MNMST is expensive (Additional file 1: Section 1.5). The comparative comparison of various algorithms on the time and space complexity is performed with different spatial transcriptomics data, where Giotto and BayesSpace are slower than stLearn and MNMST, and SCANPY, SpaGCN, and DeepST are faster than others (Additional file 1: Section 1.6 and Fig. S4). Specifically, running time of MNMST for DLPFC data is 4.3 ± 2.5 minutes, while that is 10.6 ± 2.9 (Giotto), 6.3 ± 2.2 (BayesSpace), 3.6 ± 0.4 (stLearn), 2.2 ± 0.4 (DeepST), 1.4 ± 0.3 (SpaGCN), 0.6 ± 0.1 (SEDR), and 0.1 ± 0.0 (SCANPY) respectively (Additional file 1: Fig. S4 A1). And, these results further demonstrate that MNMST achieves the best performance by reaching a good tradeoff between space and running time, proving its superiority for spatial domain identification. Furthermore, MNMST can also be accelerated with hardware, extending its applicability to large data (Additional file 1: Section 1.6 and Fig. S4).

MNMST constructs the cell affinity network with self-representation learning (SRL) and sparsity constraints, whereas baselines directly utilize cell co-expression networks by calculating the similarity of cells [16, 25]. Then, we check how performance of MNMST changes by replacing the learned cell affinity network with five cell co-expression networks (Spearman, Pearson, KNN, Cosine, and Euclidean distance). Figure 3A depicts ARIs of MNMST with various cell networks on DLPFC data (SRL: 0.553 ± 0.073 vs KNN: 0.446 ± 0.069, *p *= 3.5E−4, Student’s *t*-test). SRL achieves the highest ARI value than others, which indicates that MNMST is more precise on characterizing the structure of spatial domains than current baselines. Furthermore, various cell co-expression networks result in distinct hierarchical structure for the same slice (151675), where ARI of Pearson is 0.611 and that of KNN is 0.435 for slice 151675. Therefore, spatial domains cannot be fully characterized with cell co-expression network. Then, we investigate the difference between cell networks learned by MNMST and those used by baselines by checking the degree distribution of networks with cumulative distribution function (CDF), which calculates probability of random variable if it is less or equal a threshold. Figure 3B demonstrates that the learned networks significantly differ from these used by baselines (p<2.2E-16, Kolmogorov-Smirnov test), accounting for why the proposed algorithm is promising for spatial domain identification.

Since MNMST jointly learns features of cells, we compare the ARI of MNMST with/without joint learning (Fig. 3C left panel), where joint learning strategy is promising for learning features of cells (joint = 0.553 ± 0.073 vs non-joint = 0.450 ± 0.078, *p*-value = 2.4E−3, Student’s *t*-test). It shows that joint learning of cell spatial and expression networks enhances quality of features of cells because it implicitly reconciles spatial and expression information. Since MNMST utilizes high-order topology structure of cell spatial network for spatial domains (details in the “Methods” section), we also compare the ARI of MNMST with 1st-order and high-order structures of networks. Right panel of Fig. 3C shows that the high-order structure of networks is significantly more precise than 1st-order for spatial domain identification (high-order: 0.553 ± 0.073 vs 1st-order: 0.488 ± 0.077, *p *= 1.7E−3, Student’s* t*-test). The reason is that the relations among cells in spatial transcriptomics data are non-linear, which cannot be fully characterized with the low-order structure.

We systematically evaluate the hyperparameters of MNMST with various spatial transcriptomics data. How dimensional of cell features effect performance of MNMST is shown in Fig. 3D, where the number of dimensions in [100, 200] is a good choice because performance of MNMST is stable (Student’s *t*-test). Furthermore, parameter effect analysis not only demonstrate that MNMST is quite stable, but also suggest the values of parameters for typical spatial transcriptomics data (Additional file 1: Fig. S5 and Section 1.7).

### Multi-layer network provides an effective and efficient strategy for integrative analysis of spatial transcriptomic data

Extensive applications of spatial sequencing technologies generate a great mount of spatially omics data, which poses a great challenge on their integrative analysis. Given that current algorithms fail to integrate spatially omics data, it is interesting to design algorithms for integrating multiple datasets from various technologies. MNMST provides a novel strategy for integrative analysis of spatial transcriptomics data from three perspectives, i.e., joining multiple sections of whole tissues, integrating multiple slides of specific tissues, and integrating datasets from various batches.

Large tissues require multiple slides to cover different sections, and we employ PASTE [34] to horizontally align the tissue slices for ensuring their spatial adjacency. Then, MNMST directly performs spatial domain identification on aligned slices. Mouse brain Visium datasets of 10 \(\times\) Genomics (https://www.10xgenomics.com/), including anterior and posterior brain slices, are selected to validate performance of algorithms, where aligned slices are shown in Fig. 3E (left top panel), and zoomed regions correspond to *CA* (cornu ammonis) and *DG* (dentate gyrus). The annotation of domains for aligned slices is from the Allen Mouse Brain Atlas (right top panel of Fig. 3E). Visualization of spatial domains obtained by DeepST and MNMST is shown in Fig. 3E (bottom panels). MNMST significantly outperforms baselines since it precisely discriminates *CA* and *DG* regions, where these domains are mixed (Additional file 1: Fig. S6). It demonstrates that MNMST preserves the structure of spatial domains across various slices, which implies that the multi-layer network model is promising for integrating spatially omics data.

Visualization of topological structure of the affinity graph of cells is shown in Fig. 3F (left top panel), where cells with the same color form a domain, to determine the reason behind the superiority of MNMST over baselines. Obviously, *CA* and *DG* correspond to two clusters in the affinity graph learned by MNMST, indicating that the topological structure of cells facilitates the identification of spatial domains. Specifically, degrees of cells in *CA* are significantly higher than those of *DG* (top middle panel of Fig. 3F, *CA*: 2.7 ± 0.4 vs *DG*: 2.4 ± 0.1, *p *= 3.4E−5, Student’s *t*-test), as well as in closeness (top right panel of Fig. 3F, *p *= 8.6E−5, Student’s *t*-test). Therefore, MNMST maintains consistence of the separation boundary of spatial domains across various slices in terms of topological structure of affinity graph. Moreover, we investigate whether MNMST can preserve the structure of domains split by various slices (surrounded by white dashed line in the right bottom panel of Fig. 3F). Visualization of the topological structure of affinity graph is shown in Fig. 3F (left bottom panel), where cells are mixed given that they belong to the same spatial domain. Topological structure of cells in the anterior and posterior slice is subtle (bottom panel of Fig. 3F). These results prove that MNMST is efficient for horizontally integrating spatially omics data.

Then, we hypothesize that vertically adjacent slices are very similar, and integrate them to facilitate spatial domain identification. The proposed multi-layer network model naturally provides an effective strategy for integration of multiple vertically adjacent slices. Specifically, MNMST first stacks multiple slices into one as input, and then learns the affinity graph from the stacked slice to promote identification of spatial domains. Experiments demonstrate that integrative analysis enhances performance of spatial domain identification on DLPFC data (151673, 151674, 151675 and 151676, Additional file 1: Fig. S7 A). In details, ARI of spatial domains after integration is 0.570 for slice 151673, whereas it is 0.554 before integration, which shows that MNMST precisely captures complementary information from adjacent slices to promote spatial domain identification.

Moreover, we also investigate capability of MNMST for vertical integration with MERFISH data, where DeepST and MNMST achieve the comparable performance, demonstrating that multi-layer network model is insensitive to platforms (ARI: 0.493 vs 0.500, Additional file 1: Fig. S7B). These experiments for integrative analysis of spatial transcriptomics data are independently performed by MNMST. Thus, it is natural to ask whether MNMST can combine with other algorithms, such as PASTE [34], to perform integrative analysis. In other words, MNMST takes data integrated by PASTE [34] as inputs and then performs spatial domain identification. Notably, MNMST outperforms DeepST on the center slice generated by PASTE, where ARI of MNMST and DeepST is 0.597 and 0.537 respectively, demonstrating that multi-layer network model is more precise to capture intrinsic structure of spatial domains (Additional file 1: Fig. S7C). Furthermore, MNMST also achieves an excellent performance on the vertical stacking slices generated by PASTE (Additional file 1: Fig. S7D).

Then, we investigate the capability of MNMST to remove batch effect from spatial transcriptomics datasets. Up-to-date methods (DeepST, stLearn, and SEDR) address these issues with additional computational modules, while MNMST can be directly applied without modifying algorithms. Specifically, MNMST first stacks multiple slices from various batches into one slice with SCANPY [16] and then constructs cell spatial network and cell expression network for the integrated slice, where cross-slice spatial proximity of spots enhances consistence of features of adjacent spots, thereby removing batch effect at some extent. With four DLPFC slices (1516173-151676), MNMST achieves the best performance, where the ARI of MNMST is 0.609 (ARI of DeepST: 0.448, ARI of stLearn: 0.248, and ARI of SEDR: 0.421, Additional file 1: Fig. S8A). Interestingly, integration of multiple slices enhances the performance of spatial domain identification, which shows that MNMST can effectively remove batch effect by retaining biological content. Except for DLPFC data, mouse slices with breast cancer are aligned vertically using PASTE, where there is obvious batch effect (Additional file 1: Fig. S8B). It is observed that spatial domains in vertically adjacent slices are highly consistent after integration, whereas these domains are mixed without integration (Additional file 1: Fig. S8). These results demonstrate that MNMST effectively leverages spatial information to correct batch effects in spatial transcriptomics data.

### Multi-layer network model precisely dissects cancer and non-cancer spatial domains

Spatial transcriptomics technologies are successfully applied to cancers, and the generalization power of the proposed multi-layer network model with cancer data can be naturally investigated. The public spatial transcriptomics data of human breast cancer generated by 10 \(\times\) Visium, consisting of 3798 spots and 36,601 genes, is selected to validate the capability of identifying cancer-related spatial domains (Fig. 4A). The dataset is manually annotated by pathologists based on H&E image and the spatial expression of reported breast cancer marker genes, which consists of 20 regions and 4 main morphotypes, i.e., ductal carcinoma in situ/lobular carcinoma in situ (DCIS/LCIS), healthy tissue (Healthy), invasive ductal carcinoma (IDC), and tumor edge (Fig. 4A). Spatial domains identified by MNMST are highly consistent with the manual annotations, with an ARI of 0.662 (Fig. 4B), whereas domains obtained by baselines exhibit less regional continuity and more outliers (Additional file 1: Fig. S9). These finding implies that multi-layer network model is also promising for characterizing and extracting cancer spatial domains.

However, state-of-the-art algorithms fail to precisely identify tumor edge regions because transition from healthy to tumor is difficult to characterize and capture, which poses a great challenge in learning discriminative features of cells. The proposed multi-layer network model captures differences among various domains because indirect topological relations among cells are exploited. Specifically, the topological structure of sub-graph induced by cells in healthy, tumor edge, and tumor domains is visualized in Fig. 4C, where spatial domains correspond to clusters in the learned affinity graph. Specifically, distributions of degrees of cells in these three domains differ greatly (Fig. 4C left, Healthy: 33.2 ± 11.7 vs Tumor: 42.0 ± 17.6, *p *= 4.1E−16, Tumor edge: 28.7 ± 10.0 vs Tumor: 42.0 ± 17.6, *p *= 6.0E−18, Student’s *t*-test), as well as in closeness (Fig. 4C right). These results further demonstrate that the multi-layer network model is also effective to characterize cancer spatial domains.

Heterogeneity of tumors is critical for the diagnosis and therapy of cancers [35]. Therefore, we hypothesize that tumor heterogeneity can also be reflected by cancer spatial domains. The Pearson correlation coefficient among domains is shown in Fig. 4D, where domains are clearly divided into two groups, i.e., tumor (domain 0, 6), and non-tumor (the rest ones). This result indicates that spatial domains can serve as bio-markers to characterize tumor heterogeneity. Moreover, we perform differentially expressed analysis and gene ontology enrichment analysis to identify differentially expressed genes (DEGs) by comparing transcriptional differences among spatial domain 0 (Healthy), spatial domain 1 (IDC), and spatial domain 2 (DCIS/LCIS) (Fig. 4F, details in the “Methods” section).

To obtain intra-tumoral transcriptional differences of genes, we detect 396 DEGs (\(\vert\)log fold change\(\vert\) \(\ge\) 2, adjusted *p *< 5.0E−2, Wilcoxon rank-sum test for significance corrected by Benjamini-Hochberg test) between domains 2 and 1 (Fig. 4E–G, Additional file 1: Fig. S10 and Additional file 2). Functions enriched by DEGs are listed in Fig. 4H, where upregulated genes are involved in immune-related pathways, and downregulated ones are associated with signal pathways of fibroblast proliferation (*p *= 4.2E5, Hyper-geometric test). *CXCL14* encodes secreted proteins involved in immuno-regulatory and inflammatory processes, which promotes tumor growth in breast cancer [36], and *COX6C* serves as a bio-marker for cancer therapy [37]. Interestingly, these genes are captured with spatial domains as shown in Fig. 4I, where *CXCL14* is upregulated in IDC, and *COX6C* is depressed in IDC. Furthermore, *APOE*, *APOC1*, *C1QA*, *C1QB*, and *NUPR1* are indicative markers of variations in the infiltration of tumor-associated macrophages (TAM) [38, 39], where TAM infiltration is known to be associated with unfavorable survival outcomes in solid tumors. This association is owing to its role in promoting tumor angiogenesis, which induces tumor migration, invasion, and metastasis [40,41,42].

Then, we investigate whether spatial domain DEGs are associated with survival time of patients by employing Kaplan-Meier survival analysis with gene expression profiles and clinical data of breast cancer in The Cancer Genome Atlas (TCGA). We find that 31 of 396 DEGs separate patients into high- and low-risk groups with significant survival time (Additional file 1: Fig. S12). For example, *MORF4L2*, *PKG1*, *VPS28*, and *KLHDC7B* predict the survival time of patients, indicating that cancer domains also facilitate the identification of bio-markers for breast cancer. For example, it is proven that *KLHDC7B* regulates interferon signaling pathway, which is critical for breast cancer tumorigenesis [43].

Analogously, we perform differential expression analysis by comparing domain 1 vs 0 and identify some interesting DEGs that are highly associated with breast cancer (Additional file 1: Fig. S10, Additional file 2,3). For example, *PRDX1* participates in the signaling pathway of fibroblast proliferation, which serves as a bio-marker to characterize progression and metastasis of human breast cancer [44]. Additionally, upregulated *B2M* is also associated with a poorer prognosis [45], and *NECTIN2* is a potential target of antibody therapy for breast cancer [46]. Furthermore, gene analysis between domains 2 and 0 indicate that *NUPR1* is upregulated, which is linked to chemotherapy resistance of breast cancer (Additional file 1: Fig. S10 and S12, Additional file 2) [47]. These results demonstrate that multi-layer network model also effectively identifies cancer spatial domains, which may shed light on revealing underlying mechanisms of cancers.

### Multi-layer network model is applicable for spatial omics data with various platforms

Except for the 10 \(\times\) Genomics Visium platform, we further determine the applicability of the proposed multi-layer network model with spatial transcriptomics data from three additional platforms, i.e., imaging-based molecular data (mFISH [48, 49]) and high-resolution spatial transcriptomics data (Stereo-seq [12] and Slide-SeqV2 [50]).

The imaging-based molecular data consist of lattice- and non-lattice-shaped structure, and the lattice-shaped STARmap data for mouse visual cortex is selected, which contains 1207 cells, 1020 genes in each cells, and 7 layers as shown in Fig. 5A. MNMST outperforms baselines with an ARI of 0.624 (Fig. 5A right, Additional file 1: Fig. S14A), whereas that of DeepST is 0.550. Notably, MNMST successfully identifies *L1* and *L2/3* layers, whilst baselines fail to discriminate these layers (surrounded by a dashed square in Fig. 5A). We further investigate the difference of cells in *L1* and *L2/3* layers in terms of topological structure of networks, where distributions of degrees of cells significantly differ (degree L1: 30.9 ± 8.6 vs L2/3: 34.1 ± 10.9, *p *= 3.6E−5, Student’s *t*-test), so does closeness. These results demonstrate that MNMST can also precisely model spatial domains in lattice-shaped STARmap data.

Given that some out-of-date algorithms, such as BaysecSpace [18], are deliberately designed for lattice-shaped spatial datasets, we evaluate performance of MNMST with non-lattice-shaped spatial transcriptomics data generated by osmFISH [49], where spatial domains are labeled with different colors (Additional file 1: Fig. S14B). MNMST substantially outperforms baselines for non-linear and non-convex spatial domains with ARI = 0.600, whereas that of DeepST 0.490 (Additional file 1: Fig. S14C). Specifically, MNMST effectively distinguishes layer 2–3 lateral and medial regions that cannot be discriminated by baselines. The visualization of topological structure of different spatial domains demonstrate that these regions are precisely modeled by the cell affinity graph learned by MNMST, and cells within different domains significantly differ in terms of topological indexes (Additional file 1: Fig. S14D). For example, significant difference of cell degree between *Lateral* and *Medial* regions is observed (*Lateral*: 2.0 ± 0.7 vs *Medial*: 2.9 ± 1.0, *p *= 2.5E−35, Student’s *t*-test). These results prove that the multi-layer network model is also effective for identifying spatial domains in imaging-based molecular platform.

Then, we apply MNMST to the Stereo-seq data from mouse olfactory bulb tissues [12], which fulfills subcellular spatial resolution using DNA nano-ball patterned array chips. It contains the rostral migratory stream (RMS), granule cell layer (GCL), internal plexiform layer (IPL), mitral cell layer (MCL), external plexiform layer (EPL), and olfactory nerve layer (ONL) (Fig. 5C left panel). MNMST precisely identifies spatial domains, where the laminar organization is consistent with annotation of layers (Fig. 5D right panel). In particular, MNMST accurately recognizes the narrow tissue structure MCL, which is validated with the expression of mitral cell marker *Scg2*. Moreover, we identify bio-marker genes for other domains, including *Pcp4* for IPL, *Olfm1* for EPL, and *Cck* for GL, respectively (Additional file 1: Fig. S14C).

Finally, we apply MNMST to a mouse olfactory bulb slice profiled by Slide-SeqV2 [50], where the annotation is from the Allen Reference Atlas [33] (Fig. 5D left panel). Spatial domains obtained by MNMST are highly consistent with annotations (Fig. 5D right panel). Specifically, MNMST accurately identifies the accessory olfactory bulb (AOB), and granular layer of the accessory olfactory bulb (AOBgr). We examine the expression levels of bio-marker genes for each domain to validate the identified spatial domains, and we find that these domains identified by MNMST are well supported by the spatial variable genes (Fig. 5E). For example, *Pcp4* is upregulated in MCL, which is consistent with previous biological experiments [51]. Collectively, these results demonstrate that the multi-layer network model precisely identifies domains in spatial domains regardless of platforms, showing the superiority of network-based model.

## Discussion

In this study, we propose a novel multi-layer network model (MNMST) to integrate gene expression and spatial information of spatial transcriptomics data for spatial domain identification. MNMST is extensively evaluated on various data from various species and tissues. Extensive experimental results demonstrate that MNMST not only precisely identifies spatial domains (Fig. 2), but also is applicable to spatial transcriptomics data generated from various platforms, such as 10 \(\times\) Genomics, mFISH, Stereo-seq, and Slide-SeqV2 (Fig. 5). Furthermore, MNMST precisely dissects cancer-related domains from cancer spatial transcriptomics data (Fig. 4).

The proposed multi-layer network model constructs the cell spatial and expression network and then performs feature learning and spatial domain identification by exploiting the structure of networks. MNMST serves as a flexible framework for spatial transcriptomics data, which can be easily extended for particular situations. First, although MNMST automatically learns cell multi-layer networks from spatial transcriptomics data, it also takes available networks as inputs. Therefore, MNMST offers users an opportunity to construct cell networks by incorporating expert-knowledge with particular backgrounds. Second, MNMST simultaneously addresses spatial transcriptomics data regardless of platforms such that users can directly apply it to available data without modifying the structure of algorithms. Third, MNMST learns an affinity graph of cells, where spatial domains correspond to clusters of the graph. It therefore improves the interpretability of spatial domains, facilitating the down-stream analysis of spatial transcriptomics data. Finally, although MNMST mainly focuses on the integration of gene expression profiles and spatial information to identify spatial domains, it can be easily extended to integrate supplementary information, such as morphological information, by constructing additional networks.

Rapid developments of spatial transcriptomics technologies can measure large number of cells with high spatial resolutions at an unprecedent speed, which consequently results in the explosion of data. However, cell multi-layer networks constructed by MNMST is proportional to the number of cells, which poses a great challenge in accelerating MNMST for large- and super large-scale cell networks with millions of cells, hampering its application to very large data. Therefore, it is one of important issues to accelerate the proposed multi-layer network model. We already show that hardware acceleration is a good choice (Additional file 1: Section 1.6 and Fig. S4C). There are multiple possible strategies to address this issue, i.e., parallel computation or distributed learning systems, as well as approximating computation. Moreover, integrating genomics and spatial omics to further enhance the resolutions and annotation of spatial domains is also interesting.

## Conclusion

In conclusion, we design a promising and novel approach MNMST that accurately identifies spatial domains with the network-based model,. This method facilitates the identification of tissue organization and enables the discovery of corresponding gene markers, providing an effective and efficient model to understand complex biological systems in a spatial context. In contrast to existing approaches, MNMST demonstrates its superiority for spatial domain identification and allows for the integrative analysis of spatial transcriptomics across multiple tissue sections. We demonstrate the benefits of MNMST through comprehensive experiments with various spatial transcriptomics data generated by different platforms, providing a novel multi-layer network model for spatial transcriptomics data.

## Methods

### Data preprocessing

For all datasets, spots (cells) outside the primary tissue regions are removed. By using SCANPY package [16], the raw expression data is filtered and normalized with log-transformed according to library sizes. By following Seurat [52], genes that are expressed in less than 10 spots are also filtered. Only top 3000 genes remain for downstream analysis according to variance of gene expression, and expression of each spot is enhanced with its adjacent ones with BANKSY [24]. Principal component analysis (PCA) is employed for dimension reduction on the augmented gene expression profiles to obtain the low-dimensional cell expression matrix \(X\in R^{n\times \iota }\), where *n* and \(\iota\) is the number of cells and dimensions respectively (usually \(\iota\)=50 [20, 25]).

### Construction of cell multi-layer networks

MNMST constructs a graph for each of spatial and expression information in spatial transcriptomics data, resulting in the multi-layer network \(\mathcal {G}=\{G^{[s]}, G^{[e]}\}\), where \(G^{[s]}\) is the cell spatial network, and \(G^{[e]}\) corresponds to the cell expression network. The adjacent matrix of a cell network *G* is denoted by \(W\in R^{n\times n}\), where element \(w_{ij}\) denotes weight on edge connecting the *i*-th and *j*-th cell. Thus, constructing a network is equivalent to obtaining its adjacent matrix.

On the construction of cell spatial network \(G^{[s]}\), MNMST first constructs a initial cell spatial network \(W^{[s]}\) with using *K*-nearest neighbor algorithm (by following Squidpy [53], the number of neighbors *k* is 6 for 10 \(\times\) Genomics data, 8 for Slide-seq and Stereo-seq data, and 15 for imaging-based platforms), where weight \(w_{ij}^{[s]}\) for the *i*-th and *j*-th cell is inverse proportional to the spatial Euclidean distance between them (denoted by \(r_{ij}\)), i.e., \(w_{ij}^{[s]}=(r_{ij})^{-1}\). Recently, evidence demonstrates that indirect relations, such as multiple hops among cells, is more precise characterize topological structure of network [54]. MNMST takes pointwise mutual information (PMI) matrix [55] as adjacent matrix of cell spatial network \(G^{[s]}\) (denoted by \(M^{[s]}\)). And, element \(m_{ij}^{[s]}\) of \(M^{[s]}\) is defined as

where \(d_l\) is the degree of l-th cell in network \(W^{[s]}\), and \(\kappa\) is the number of non-negative samples (default set to 1). In this study, we refer PMI matrix \(M^{[s]}\) as the high-order topological structure of network *G*.

On the construction cell expression network \(G^{[e]}\), the cell expression matrix *X* is first augmented with the initial cell spatial network \(W^{[s]}\) according to Ref. [24]. Then, MNMST utilizes sparse self-representation learning (SRL) [56,57,58] with local preservation constraint [59] to automatically learn adjacent matrix \(W^{[e]}\), which is under the assumption that each cell can be well represented with its neighbors. And, the objective function for constructing cell expression network \(G^{[e]}\) is formulated as a minimization problem as

where \(L^{[e]}=D^{[e]}-W^{[e]}\) is the Laplacian matrix, *Tr*(*A*) denotes the trace of matrix *A* (i.e., \(Tr(A)=\sum _{i}a_{ii}\)), \(\Vert \cdot \Vert ^2\) and \(\Vert \cdot \Vert _1\) are \(l_2-\) and \(l_1\)-norm respectively, and \(\alpha\) and \(\beta\) are parameters. The constraint \(diag(W^{[e]})=0\) avoids trivial solutions (formulation and optimization of Eq. 2 are shown in Additional file 1: Section 1.1 and 1.2).

### Spatial and transcriptional feature learning by factorizing cell multi-layer networks

Given cell multi-layer networks, the most intuitive strategy is to independently learn cell features for each graph, which are adopted by current algorithms because of simplicity. However, this strategy is criticized for the low quality of cell features considering that it neglects relations between cell spatial and expression network. In our previous study [60], we demonstrate that jointly factorizing multi-layer networks is promising in capturing the relations of various graph, which significantly enhances the quality of features. Herein, MNMST utilizes a joint model with non-negative matrix factorization [61] to decompose the cell multi-layer network as

where matrix *B* denotes the shared feature of cells.

Similar to Eq. 2, MNMST automatically learns an affinity graph based on matrix *B* with the denoising model and self-representation learning, which is formulated as

where \(B'\) denotes the transpose of *B*, *E* is the error matrix, \(\Vert Z\Vert _*\) is the nuclear norm of *Z* [62], and \(\Vert E\Vert _{2,1}\) is the \(l_{2,1}\)-norm, respectively (mathematical model for Eq. 4 in Additional file 1: Section 1.1).

Independence of feature learning and affinity graph construction fails to characterize patterns in spatial omics data, resulting in unsatisfactory spatial domain identification performance entirely, i.e., MNMST without joint learning. Previous studies [15] demonstrate that joint learning is more precise to characterize structure of networks, and [63] verifies that joint learning is also applicable to spatial omics data. Inspired by previous studies, MNMST joints feature learning and affinity graph construction by combining Eqs. 3 and 4; objective function of feature learning is formulated as

where \(\gamma\) and \(\lambda\) are parameters for tuning. The Alternating Direction Method of Multipliers (ADMM) method [64] is employed to optimize Eq. 5, and update rules are deduced (Additional file 1: Section 1.3).

### Clustering and visualization

Based on the learned affinity graph, MNMST employs Leiden algorithm [65] to obtain graph partitioning, where each cluster corresponds to a spatial domain. If the number of spatial domains is known in advance, MNMST performs grid search with gap of 0.01 until the number of clusters is reached. Instability of nonnegative matrix factorization is adopted to obtain the number of domains when it is unknown [60]. To visualize distributions of cells in tissues, the uniform manifold approximation and projection (UMAP) algorithm is selected [66].

### Clustering criteria

When annotations of spatial domains are missing, two extensively adopted clustering criteria, including Silhouette Coefficient (SC) and Davies-Bouldin (DB) scores, are selected to validate the performance of clustering in terms of computation. Specifically, SC takes into account compactness within and separation across clusters as \((b-a) / max(a, b)\), where *a* is the mean intra-cluster distance, and *b* is the mean nearest-cluster distance. It ranges between −1 and 1, where a higher score refers to more coherent clusters. SC=0 means that the sample is on or close to the boundary of neighboring clusters, whereas negative values denote potentially wrong clusters. DB score is the average ratio of within-cluster distances to between-cluster distances, favoring farther apart and less dispersed clusters with low values.

### Spatial trajectory inference

PAGA [32] infers trajectories of cell types for single-cell transcriptomics data, which generates graph-like maps of cells by preserving continuous and disconnected structures at multiple resolutions. Specifically, it outputs a graph, where nodes are clusters and edges are connectivity (or similarity) of clusters. Here, we adopt it to infer relations of spatial domains by inputting affinity graph of cells as well as labels of domains.

### Identification and functional analysis of differentially expressed genes

MNMST performs differential expression analysis of genes for each spatial domain by using Wilcoxon rank-sum test implemented in SCANPY package [16]. Genes are expressed in 80% cells/spots in each domain and are with the absolute value of log fold change \(\ge\) 2, and an adjusted FDR \(\le\) 0.05 is selected as differentially expressed genes (DEGs). The overlapping genes between DEGs and highly variable genes (HVG) obtained by spatialDE [27] are used for down-stream analysis. Gene ontology enrichment analysis is performed with clusterProfiler [67].

### Benchmarking

Seven state-of-the-art methods, namely SCANPY [16], Giotto [29], stLearn [28], SEDR [22], BayesSpace [18], SpaGCN [20], and DeepST [25], are selected as baselines, where SCANPY is non-spatial, and others are spatial clustering algorithms. All these algorithms are executed with the suggested values of parameters to achieve the best performance for fair comparison.

For datasets with known annotations of spatial domains, given ground truth \(P^{*}\) and predicted domains *P*, the adjusted rank index (ARI) is defined as [31]

where *n* is the number of cells, \(n_{ij}\) is the number of cells of class label \(C^{*} \in P^{*}\) assigned to cluster \(C_i\) in partition *P*, and \({n_i}\) / \({n_j}\) is the number of cells in cluster \(C_i\) / \(C_j\) of partition *P*.

## Availability of data and materials

The code for the MNMST algorithm is implemented in Python and detailed tutorials are freely available at GitHub [68] as well as Zenodo [69]. The source code is released under the MIT License. All datasets used in this paper are published datasets available for download. Human DLPFCs is accessible within the SpatialLIBD [70] at http://spatial.libd.org/spatialLIBD. Human breast cancer and mouse brain tissue sections datasets are collected from the 10 \(\times\) Genomics website at https://www.10xgenomics.com/resources/datasets. The mouse visual cortex STARmap data [71] is accessible at https://www.wangxiaolab.org/data-portal-1. The mouse brain cortex osmFISH data [72] is accessible at http://linnarssonlab.org/osmFISH. Slide-seqV2 dataset [73] is available at the Broad Institute Single Cell Portal at https://singlecell.broadinstitute.org/single_cell/study/SCP815. The processed Stereo-seq data [74] from mouse olfactory bulb tissue is accessible at https://github.com/JinmiaoChenLab/SEDR_analyses. The MERFISH data [75] collected from the mouse hypothalamic preoptic region is accessible at https://zhuang.harvard.edu/merfish.

## References

Rao A, Barkley D, França GS, Yanai I. Exploring tissue architecture using spatial transcriptomics. Nature. 2021;596:211–20.

Zhang AW, O’Flanagan C, Chavez EA, Lim JL, Ceglia N, McPherson A, et al. Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling. Nat Methods. 2019;16:1007–15.

Andrews N, Serviss JT, Geyer N, Andersson AB, Dzwonkowska E, Šutevski I, et al. An unsupervised method for physical cell interaction profiling of complex tissues. Nat Methods. 2021;18:912–20.

Li B, Zhang W, Guo C, Xu H, Li L, Fang M, et al. Benchmarking spatial and single-cell transcriptomics integration methods for transcript distribution prediction and cell type deconvolution. Nat Methods. 2022;19:662–70.

Asp M, Bergenstråhle J, Lundeberg J. Spatially resolved transcriptomes-next generation tools for tissue exploration. Bioessays. 2020;42:e1900221.

Williams CG, Lee HJ, Asatsuma T, Vento-Tormo R, Haque A. An introduction to spatial transcriptomics for biomedical research. Genome Med. 2022;14:1–18.

Levsky JM, Singer RH. Fluorescence in situ hybridization: past, present and future. J Cell Sci. 2003;116:2833–8.

Moffitt JR, Zhuang X. RNA imaging with multiplexed error-robust fluorescence in situ hybridization (MERFISH). In: Methods in Enzymology. vol. 572. Elsevier; 2016. p. 1–49.

Moffitt JR, Hao J, Wang G, Chen KH, Babcock HP, Zhuang X. High-throughput single-cell gene-expression profiling with multiplexed error-robust fluorescence in situ hybridization. Proc Natl Acad Sci. 2016;113:11046–51.

Moffitt JR, Bambah-Mukku D, Eichhorn SW, Vaughn E, Shekhar K, Perez JD, et al. Molecular, spatial, and functional single-cell profiling of the hypothalamic preoptic region. Science. 2018;362(6416):eaau5324.

Rodriques SG, Stickels RR, Goeva A, Martin CA, Murray E, Vanderburg CR, et al. Slide-seq: a scalable technology for measuring genome-wide expression at high spatial resolution. Science. 2019;363:1463–7.

Wei X, Fu S, Li H, Liu Y, Wang S, Feng W, et al. Single-cell Stereo-seq reveals induced progenitor cells involved in axolotl brain regeneration. Science. 2022;377:eabp9444.

Ståhl PL, Salmén F, Vickovic S, Lundmark A, Navarro JF, Magnusson J, et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science. 2016;353:78–82.

Zhang M, Eichhorn SW, Zingg B, Yao Z, Cotter K, Zeng H, et al. Spatially resolved cell atlas of the mouse primary motor cortex by MERFISH. Nature. 2021;598:137–43.

Wu W, Ma X. Joint learning dimension reduction and clustering of single-cell RNA-sequencing data. Bioinformatics. 2020;36:3825–32.

Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19:1–5.

Blondel VD, Guillaume JL, Lambiotte R, Lefebvre E. Fast unfolding of communities in large networks. J Stat Mech Theory Exp. 2008;2008:P10008.

Zhao E, Stone MR, Ren X, Guenthoer J, Smythe KS, Pulliam T, et al. Spatial transcriptomics at subspot resolution with BayesSpace. Nat Biotechnol. 2021;39:1375–84.

Shang L, Zhou X. Spatially aware dimension reduction for spatial transcriptomics. Nat Commun. 2022;13:7203.

Hu J, Li X, Coleman K, Schroeder A, Ma N, Irwin DJ, et al. SpaGCN: Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network. Nat Methods. 2021;18:1342–51.

Dong K, Zhang S. Deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder. Nat Commun. 2022;13:1–12.

Xu H, Fu H, Long Y, Ang KS, Sethi R, Chong K, et al. Unsupervised spatially embedded deep representation of spatial transcriptomics. Genome Med. 2024;16(1):12.

Li J, Chen S, Pan X, Yuan Y, Shen HB. Cell clustering for spatial transcriptomics data with graph neural networks. Nat Comput Sci. 2022;2:399–408.

Singhal V, Chou N, Lee J, Yue Y, Liu J, Chock WK, et al. BANKSY unifies cell typing and tissue domain segmentation for scalable spatial omics data analysis. Nat Genet. 2024:1–11.

Xu C, Jin X, Wei S, Wang P, Luo M, Xu Z, et al. DeepST: identifying spatial domains in spatial transcriptomics by deep learning. Nucleic Acids Res. 2022;50:e131.

Sun S, Zhu J, Zhou X. Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies. Nat Methods. 2020;17:193–200.

Svensson V, Teichmann SA, Stegle O. SpatialDE: identification of spatially variable genes. Nat Methods. 2018;15:343–6.

Pham D, Tan X, Balderson B, Xu J, Grice LF, Yoon S, et al. Robust mapping of spatiotemporal trajectories and cell–cell interactions in healthy and diseased tissues. Nat Commun. 2023;14(1):7739.

Dries R, Zhu Q, Dong R, Eng CHL, Li H, Liu K, et al. Giotto: a toolbox for integrative analysis and visualization of spatial expression data. Genome Biol. 2021;22:1–31.

Maynard KR, Collado-Torres L, Weber LM, Uytingco C, Barry BK, Williams SR, et al. Transcriptome-scale spatial gene expression in the human dorsolateral prefrontal cortex. Nat Neurosci. 2021;24:425–36.

Hubert L, Arabie P. Comparing partitions J Classif. 1985;2:193–218.

Wolf FA, Hamey FK, Plass M, Solana J, Dahlin JS, Göttgens B, et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. 2019;20:1–9.

Sunkin SM, Ng L, Lau C, Dolbeare T, Gilbert TL, Thompson CL, et al. Allen Brain Atlas: an integrated spatio-temporal portal for exploring the central nervous system. Nucleic Acids Res. 2012;41:D996–1008.

Zeira R, Land M, Strzalkowski A, Raphael BJ. Alignment and integration of spatial transcriptomics data. Nat Methods. 2022;19:567–75.

Ge R, Wang Z, Cheng L. Tumor microenvironment heterogeneity an important mediator of prostate cancer progression and therapeutic resistance. NPJ Precis Onco. 2022;6:31.

Sjöberg E, Meyrath M, Milde L, Herrera M, Lövrot J, Hägerstrand D, et al. A novel ACKR2-dependent role of fibroblast-derived CXCL14 in epithelial-to-mesenchymal transition and metastasis of breast cancer. Clin Cancer Res. 2019;25:3702–17.

Jang SC, Crescitelli R, Cvjetkovic A, Belgrano V, Olofsson Bagge R, Sundfeldt K, et al. Mitochondrial protein enriched extracellular vesicles discovered in human melanoma tissues can be detected in patient plasma. J Extracell Vesicles. 2019;8:1635420.

Wang YF, Yu L, Hu ZL, Fang YF, Shen YY, Song MF, et al. Regulation of CCL2 by EZH2 affects tumor-associated macrophages polarization and infiltration in breast cancer. Cell Death Dis. 2022;13:748.

Nalio Ramos R, Missolo-Koussou Y, Gerber-Ferder Y, Bromley CP, Bugatti M, Núñez NG, et al. Tissue-resident FOLR2+ macrophages associate with CD8+ T cell infiltration in human breast cancer. Cell. 2022;185:1189–1207.e25.

Ruffell B, Coussens LM. Macrophages and therapeutic resistance in cancer. Cancer Cell. 2015;27:462–72.

Lewis CE, Pollard JW. Distinct role of macrophages in different tumor microenvironments. Cancer Res. 2006;66:605–12.

Engblom C, Pfirschke C, Pittet MJ. The role of myeloid cells in cancer therapies. Nat Rev Cancer. 2016;16:447–62.

Li MX, Wang HY, Yuan CH, Ma ZL, Jiang B, Li L, et al. KLHDC7B-DT aggravates pancreatic ductal adenocarcinoma development via inducing cross-talk between cancer cells and macrophages. Clin Sci. 2021;135:629–49.

Turner-Ivey B, Manevich Y, Schulte J, Kistner-Griffin E, Jezierska-Drutel A, Liu Y, et al. Role for Prdx1 as a specific sensor in redox-regulated senescence in breast cancer. Oncogene. 2013;32:5302–14.

Pitt JJ, Riester M, Zheng Y, Yoshimatsu TF, Sanni A, Oluwasola O, et al. Characterization of Nigerian breast cancer reveals prevalent homologous recombination deficiency and aggressive molecular features. Nat Commun. 2018;9:4181.

Oshima T, Sato S, Kato J, Ito Y, Watanabe T, Tsuji I, et al. Nectin-2 is a potential target for antibody therapy of breast and ovarian cancers. Mol Cancer. 2013;12:1–13.

Mu Y, Yan X, Li D, Zhao D, Wang L, Wang X, et al. NUPR1 maintains autolyso somal efflux by activating SNAP25 transcription in cancer cells. Autophagy. 2018;14:654–70.

Wang X, Allen WE, Wright MA, Sylwestrak EL, Samusik N, Vesuna S, et al. Three-dimensional intact-tissue sequencing of single-cell transcriptional states. Science. 2018;361:eaat5691.

Codeluppi S, Borm LE, Zeisel A, La Manno G, van Lunteren JA, Svensson CI, et al. Spatial organization of the somatosensory cortex revealed by osmFISH. Nat Methods. 2018;15:932–5.

Stickels RR, Murray E, Kumar P, Li J, Marshall JL, Di Bella DJ, et al. Highly sensitive spatial transcriptomics at near-cellular resolution with Slide-seqV2. Nat Biotechnol. 2021;39:313–9.

Kim EE, Shekhar A, Lu J, Lin X, Liu FY, Zhang J, et al. PCP4 regulates Purkinje cell excitability and cardiac rhythmicity. J Clin Investig. 2014;124:5027–36.

Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. 2015;33:495–502.

Palla G, Spitzer H, Klein M, Fischer D, Schaar AC, Kuemmerle LB, et al. Squidpy: a scalable framework for spatial omics analysis. Nat Methods. 2022;19:171–8.

Benson AR, Gleich DF, Leskovec J. Higher-order organization of complex networks. Science. 2016;353(6295):163–6.

Li Y, Sha C, Huang X, Zhang Y. Community detection in attributed graphs: an embedding approach. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 32. Palo Alto: AAAI Press; 2018.

Kolesnikov A, Zhai X, Beyer L. Revisiting self-supervised visual representation learning. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. Long Beach: IEEE Press; 2019. p. 1920–29.

Liu G, Lin Z, Yan S, Sun J, Yu Y, Ma Y. Robust recovery of subspace structures by low-rank representation. IEEE Trans Pattern Anal Mach Intell. 2012;35:171–84.

Elhamifar E, Vidal R. Sparse subspace clustering: algorithm, theory, and applications. IEEE Trans Pattern Anal Mach Intell. 2013;35(11):2765–81.

Fang X, Xu Y, Li X, Fan Z, Liu H, Chen Y. Locality and similarity preserving embedding for feature selection. Neurocomputing. 2014;128:304–15.

Ma X, Dong D, Wang Q. Communitydetection in multi-layer networks using joint nonnegative matrix factorization. IEEE Trans Knowl Data Eng. 2018;31:273–86.

Liu J, Wang C, Gao J, Han J. Multi-view clustering via joint nonnegative matrix factorization. In: Proceedings of the 2013 SIAM International Conference on Data Mining. Philadelphia; SIAM; 2013. p. 252–60.

Cai JF, Candès EJ, Shen Z. A singular value thresholding algorithm for matrix completion. SIAM J Optim. 2010;20(4):1956–82.

Liu W, Liao X, Yang Y, Lin H, Yeong J, Zhou X, et al. Joint dimension reduction and clustering analysis of single-cell RNA-seq and spatial transcriptomics data. Nucleic Acids Res. 2022;50(12):e72.

Boyd S, Parikh N, Chu E, Peleato B, Eckstein J, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found Trends Mach Learn. 2011;3:1–122.

Traag VA, Waltman L, Van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep. 2019;9:5233.

Becht E, McInnes L, Healy J, Dutertre CA, Kwok IW, Ng LG, et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2019;37:38–44.

Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2:100141.

Wang Y, Liu Z, Ma X. MNMST: Topology of cell networks leverages identification of spatial domains from spatial transcriptomics data. Github. 2024. https://github.com/xkmaxidian/MNMST. Accessed 21 May 2024.

Wang Y, Liu Z, Ma X. MNMST: Topology of cell networks leverages identification of spatial domains from spatial transcriptomics data. Zenodo. 2024. https://doi.org/10.5281/zenodo.10369465.

Maynard KR, Collado-Torres L, Weber LM, Uytingco C, Barry BK, Williams SR, et al. Transcriptome-scale spatial gene expression in the human dorsolateral prefrontal cortex. DLPFC dataset by 10\(\times\) Visium. spatiaLIBD. 2021. https://research.libd.org/spatialLIBD/. Accessed 21 May 2024.

Wang X, Allen WE, Wright MA, Sylwestrak EL, Samusik N, Vesuna S, et al. Three-dimensional intact-tissue sequencing of single-cell transcriptional states. Mouse medial prefrontal cortex dataset by STARmap. STARmap Resour. 2018. http://starmapresources.org/data. Accessed 21 May 2024.

Codeluppi S, Borm LE, Zeisel A, La Manno G, van Lunteren JA, Svensson CI, et al. Spatial organiza cortex revealed by osmFISH. Mouse visual cortex dataset by osmFISH. osmFISH Resour. 2018. http://linnarssonlab.org/osmFISH. Accessed 21 May 2024.

Stickels RR, Murray E, Kumar P, Li J, Marshall JL, Di Bella DJ, et al. Highly sensitive spatial transcriptomics at near-cellular resolution with Slide-seqV2. Mouse Olfactory bulb dataset by Slide-seq V2. Single Cell Portal. 2020. https://singlecell.broadinstitute.org/singlecell/study/SCP815. Accessed 20 Feb 2024.

Xu H, Fu H, Long Y, Ang KS, Sethi R, Chong K, et al. Unsupervised spatial embedded deep representation of spatial transcriptomics. Mouse Olfactory bulb dataset by Stereo-seq. Github. 2021. https://github.com/JinmiaoChenLab/SEDR_analyses. Accessed 21 May 2024.

Moffitt JR, Bambah-Mukku D, Eichhorn SW, Vaughn E, Shekhar K, Perez JD, et al. Molecular, spatial, and functional single-cell profiling of the hypothalamic preoptic region. Mouse hypothalamus dataset by MERFISH. DRYAD. 2018. https://doi.org/10.5061/dryad.8t8s248. Accessed 21 May 2024.

## Acknowledgements

The authors thank the members of the Ma Lab for helpful discussion and appreciate the researchers who provide us with source code for a comparison. The authors thank reviewers for their time and suggestions.

### Review history

The review history is available as Additional file 4.

### Peer review information

Ewa Szczurek and Veronique van den Berghe were the primary editors of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

## Funding

This work was supported by the Joint Funds of the National Natural Science Foundation of China (U22A20345), National Natural Science Foundation of China (62272361), and Shaanxi Natural Science Funds for Distinguished Young Scholar (Program No. 2022JC-38).

## Author information

### Authors and Affiliations

### Contributions

X.M. conceived and designed the study; X.M., Z.L., and Y.W. performed the research; X.M. and Y.W. collected and constructed the benchmark datasets and constructed the models; X.M., Z.L., and Y.W. completed the downstream analysis.

### Corresponding author

## Ethics declarations

### Ethics approval and consent to participate

No ethical approval was required for this study. All utilized public datasets were generated by other organizations that obtained ethical approval.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

## Additional information

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

### 13059_2024_3272_MOESM1_ESM.pdf

Additional file 1: Supplementary information. Supplementary notes describing the detailed derivations of the MNMST algorithm, and Supplementary figures.

### 13059_2024_3272_MOESM2_ESM.xlsx

Additional file 2: Supplementary table, differentially expressed genes (DEGs) identified between Domain 2 and Domain 1, Domain 2 and Domain 0, and Domain 1 and Domain 0 in human breast cancer dataset.

### 13059_2024_3272_MOESM3_ESM.xlsx

Additional file 3: Supplementary table, gene ontology enrichment analysis results corresponding for DEGs in Additional file 2.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

## About this article

### Cite this article

Wang, Y., Liu, Z. & Ma, X. MNMST: topology of cell networks leverages identification of spatial domains from spatial transcriptomics data.
*Genome Biol* **25**, 133 (2024). https://doi.org/10.1186/s13059-024-03272-0

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s13059-024-03272-0