A novel multislice framework for precision 3D spatial domain reconstruction and disease pathology analysis

Overview of the STMSC framework

STMSC is a deep learning data analysis framework tailored for ST data. This framework takes as input gene expression data, spatial coordinates, and histological images from ST, integrating corresponding SC data to provide a comprehensive latent representation of spatial information for subsequent biological analysis. First, STMSC performs the task of slice alignment (Fig. 1, step 1). This step begins by digitizing H&E-stained images and combining spatial coordinates with image information to achieve precise slice alignment, resulting in 3D spatial coordinates. Utilizing these coordinates, k-nearest neighbor (K-NN) graphs are constructed. Next, STMSC uses SC data to deconstruct the microenvironment, accurately capturing biological information at the cell level for each spot and creating a complete cellular landscape (Fig. 1, step 2). This information is used to further refine the neighborhood graph (Fig. 1, step 3), effectively addressing issues that may arise with traditional spatial coordinate distance measurements at complex domain boundaries. Finally, the corrected 3D neighborhood graph and ST gene expression data are input into a GAE (Fig. 1, step 4). This model learns the latent spatial representation of spots, which reflects the biological information and structural characteristics within the tissue. Once trained, these latent representations can be used for spatial domain identification, 3D tissue structure reconstruction, and various downstream analysis tasks (Fig. 1, step 5). Among the downstream analyses, trajectory inference was conducted with particular attention to biological plausibility. To this end, the human cerebral cortex data set was selected, as it follows a well-characterized linear developmental trajectory, in contrast to more structurally complex tissues such as breast cancer. Figure 1 depicts the comprehensive architectural framework of STMSC, and the Methods section provides detailed algorithmic specifics.

Figure 1.

Overview of STMSC. (1) Slice alignment: uses iterative closest point (ICP) (Arun et al. 1987) to align multiple slices, establishing the three-dimensional (3D) positions of spots and constructing a global 3D structure. (2) Microenvironmental deconstruction: trains a cell-to-spot mapping matrix M using spatial transcriptomic (ST) gene expression data (X_st) and single-cell data (X_sc) via a spatially informed contrastive learning model. In this model, the similarity of positive pairs (i.e., spatially adjacent spots) is maximized, whereas the similarity of negative pairs (i.e., spatially nonadjacent spots) is minimized. (3) Construction and correction of 3D neighborhood graph. (4) Joint modeling: STMSC jointly models multiple slices and employs a graph attention mechanism-based autoencoder to learn latent spot representations with 3D spatial information. (5) Downstream analysis: The learned latent representations serve as inputs for downstream analysis, including 3D spatial domain identification and spatial trajectory inference.

Open in new tabDownload PowerPointLink to figure

In this section, we aim to validate the efficacy of the STMSC model in accurately identifying 3D spatial domains and reconstructing complex 3D spatial structures. First, we demonstrate the model's ability to precisely delineate cortical regions and reconstruct the 3D spatial architecture of the human brain. Next, we demonstrate STMSC's ability to capture multilayered biological heterogeneity, thereby revealing subtle biological textures. Finally, we show how accurate spatial domain identification by STMSC provides deeper insights into complex cancer microenvironments. Additionally, we validate our approach using simulated data sets (Supplemental Table 1), with detailed analyses provided in Supplemental Section 2.

STMSC reconstructs the 3D spatial organization of the human brain and accurately identifies spatial domains

In this study, we utilize a data set from the dorsolateral prefrontal cortex (DLPFC) of three adult human samples (Maynard et al. 2021), each consisting of four consecutive tissue slices, generated using the 10x Genomics Visium platform (Maynard et al. 2021). These sections are meticulously annotated with different cortical layers (layers I through VI) and white matter (WM) regions. These annotations serve as a benchmark for comprehensively assessing the spatial domain recognition capabilities of the STMSC model in both single-slice and multislice data integration scenarios. The correspondence between slice names and slice numbers is provided in Supplemental Table 7.

First, the multislice alignment step, incorporating both coordinate information and histological images, effectively reconstructed 3D tissue structures (Fig. 2A), enriching the information available for 3D spatial domain identification. We then quantitatively evaluated the spatial domain recognition capabilities of STMSC using the adjusted Rand index (ARI) (Fig. 2B) and conducted a detailed comparison with existing advanced ST analysis methods, including SpaGCN (Hu et al. 2021), STAGATE (Dong and Zhang 2022), STAligner (Zhou et al. 2023), GraphST (Long et al. 2023), STitch3D (Wang et al. 2023), and SPACEL (Xu et al. 2023). In Figure 2B, the height of each bar represents the average ARI across 12 tissue slices for a given method, and each dot denotes the ARI score from an individual slice, providing a direct visualization of both overall performance and variability across slices. Notably, STMSC consistently achieves higher mean ARI values than all other methods across the 12 slices, indicating its superior performance in spatial domain identification. These gains likely result from STMSC's ability to integrate latent spatial information across tissue sections (Fig. 2D; Haghverdi et al. 2016; Wolf et al. 2019). The latent representations of all sections exhibit consistent patterns, with trajectories from the WM to layer I aligning with the neurogenesis sequence during cortical development (Gilmore and Herrup 1997). STMSC further validates the consistency between spatial domain trajectories and spatial maps, highlighting its exceptional performance in 3D spatial domain recognition.

Figure 2.

3D reconstruction of the human brain and comparison with other advanced methods. (A) 3D spatial domain identification results for DLPFC data. (B) Comparison results of clustering index ARI. (C) Mean results of ablation experiments on the STMSC. (D) UMAP visualization of STMSC's learned spot representations in the shared latent space. (E) Results of ablation experiments on the STMSC. (F) 2D visualization results of spatial domain recognition. (G) Spatial distribution of cell types on section Sample3-Slice1.

Open in new tabDownload PowerPointLink to figure

Second, we conducted a detailed analysis on four consecutive slices from Sample 1 (Fig. 2F), utilizing UMAP to visualize the spatial domain recognition results of STMSC, SpaGCN, GraphST, SPACEL, and Stitch3D. The results indicate that SpaGCN faces challenges in processing cortical boundaries, making them less clearly identifiable. GraphST shows some inconsistencies in layer integrity, leading to mixing within the same domain. SPACEL detects an additional cortical region in the upper right corner of Sample1-Slice3 and Sample1-Slice4, where such a structure is not typically found. SPACEL demonstrates better performance in identifying thinner layers, likely owing to its ability to generate smoother boundaries in structured cortical regions. Stitch3D exhibits domain mixing in spatial domain 0, which affects the accurate identification of layers. Compared with other methods, STMSC excels in accurately identifying each cortical layer (i.e., spatial domain) with minimal interlayer confusion, achieving a higher level of completeness in layer delineation than competing approaches. This outcome validates the effectiveness of STMSC's design in reconstructing 3D tissue structures, characterizing the cellular landscape, and capturing subtle biological textures within tissues, ultimately enhancing the accuracy of spatial domain identification.

Furthermore, to evaluate the reliability of microenvironmental structures identified with our framework, we analyzed sample3-slice1 as an example. Because of the lack of direct experimental evidence detailing the cellular composition at each point in this data set, we used the correlation annotations from the original study (Maynard et al. 2021) as a reference. The results (Fig. 2G) show that oligodendrocyte (Oligo) subtypes are predominantly expressed in the WM, whereas astrocyte (Astro) subtypes are preferentially expressed in layer L1. Oligodendrocyte precursor cells (OPCs), which are immature Oligos that can differentiate into mature Oligos, are located along the edges of the WM, aligning with biological expectations. Additionally, inhibitory (Inh) neurons are primarily distributed across layers L2 to L6. To further investigate, we performed a detailed visualization of excitatory neurons. According to the literature, excitatory neurons (Exc) marked by RORB and CMAHP are mainly distributed in layers L3 to L5, consistent with our observations. In conclusion, our experimental findings confirm the reliability of the microenvironmental structure module within our framework, as its analytical outcomes closely correspond to the cortical annotations presented in the original study. This alignment not only attests to the accuracy of our method but also facilitates new advancements in the study of the human cerebral cortex. By offering deeper insights into the cellular composition and spatial distribution within the cortex, our research establishes a solid foundation for future investigations in the field of neuroscience.

Finally, we conducted an ablation analysis using the DLPFC data set. We systematically evaluated the impact of each module within the STMSC framework by generating several ablated versions: baseline GAT+JA+, which removes the correction module; baseline GAT+CM+, which excludes the multislice joint analysis module; baseline GAT+, which omits both the correction mechanism and the multislice joint analysis module; and baseline, which eliminates the correction mechanism, multislice joint analysis, and GAE. The results (Fig. 2C,E) demonstrate that the average performance of STMSC across the 12 slices outperforms all ablated versions, confirming the necessity and effectiveness of each module. This indicates that adding H&E image data offers essential phenotypic information, which improves the slice alignment task and increases the effectiveness of multilayer joint analysis. The reconstruction of the cellular landscape and the insights gained into collective cellular behaviors within tissues are essential for accurate domain identification. Moreover, the introduced precorrection mechanism effectively aids in capturing subtle biological textures within tissues. In summary, the results of this experiment validate the concepts discussed earlier in this paper.

STMSC enhances domain identification by capturing multilevel biological heterogeneity, thereby revealing subtle biological textures

Accurately identifying and characterizing small spatial domains is critical in the study of biological tissue microstructures (Kalita-de Croft et al. 2021; Zheng et al. 2023). Given the highly complex architecture of the brain cortex, adjacent regions may differ in cellular composition and function. Therefore, constructing a detailed cellular landscape of the tissue is essential. STMSC addresses the challenge of missing microlevel information by employing a microenvironment deconstruction task, thereby providing critical support for domain identification. Additionally, traditional distance-based graph structuring methods may fail to define domain boundaries effectively. To address this, the STMSC model incorporates a precorrection mechanism that utilizes fine-grained cellular composition data to enhance domain delineation. The effectiveness of this mechanism and the model's capability in identifying subtle biological structures were validated using 10x Visium data from coronal brain sections of mice, annotated with gene expression maps from the Allen Brain Atlas and Mouse Brain available through Squidpy (Fig. 3A, left; Palla et al. 2022).

Figure 3.

Analysis of the mouse coronal brain data set. (A) 2D visualization of STMSC spatial domain identification versus true labels. (B) Alluvial plot showing correspondence between STMSC domains and manual annotations (ground truth). (C) 2D visualization of spatial domain identification by other methods. (D) ARI comparison of clustering on the mouse prefrontal cortex. (E) Ablation study results. (F) Dotplot of spatially differentially expressed genes.

Open in new tabDownload PowerPointLink to figure

Initially, experimental results (Fig. 3A,C) demonstrated that the spatial regions identified by the STMSC model align closely with manually annotated regions in Squidpy. Notably, within the hippocampus, STMSC precisely differentiated between two distinct structures: the pyramidal layer and the dentate gyrus. In contrast, other models like SpaGCN, GraphST, and SPACEL failed to recognize these subtle structural differences. Although the STitch3D model could roughly delineate these areas, it lacked the precision achieved by STMSC. Further analysis revealed that STMSC's ARI reached 0.66, leading the field and showcasing the model's superiority in spatial region partitioning accuracy (Fig. 3D). Figure 3B illustrates the correspondence between different STMSC spatial regions (e.g., Cortex_1, Cortex_2, Hippocampus, Thalamus_1, etc.) and the manually annotated regions. Each flow represents the mapping from the STMSC-identified region to the manual annotation labels. Different colors represent the regions, and the width of the flow indicates the degree of correspondence. The STMSC method shows good spatial domain identification ability, with high consistency in most regions compared with the manual annotations. The ablation experiment results presented in Figure 3E further demonstrate the effectiveness of the precorrection mechanism.

Subsequently, in exploring the differential expression of spatially variable genes, we confirmed that specific gene expression patterns across different brain areas are closely associated with neurodevelopment and functional disorders (Fig. 3F). For example, mutations in the Mef2c gene in a mouse model of autism affect neurodevelopment and core behavioral phenotypes, primarily evident in cortical development (Zhang and Zhao 2022; Li et al. 2024). Similarly, knockout of the Nrsn1 gene alters cortical morphology, such as thickening of the cortex and loose cell arrangement, which could impact brain functional development (Li et al. 2016; Skeide et al. 2016). Dclk1 gene expression is crucial for neuronal migration and the formation of cortical structures, and its anomalies may be linked to developmental disorders like epilepsy and autism spectrum disorders (Ding et al. 2023). Moreover, although the Sparc gene is primarily associated with obesity and metabolic functions, its role in regulating related hypothalamic functions cannot be overlooked (Chaiyasoot et al. 2023). Expression of the Spink8 gene in mouse hippocampal pyramidal cells reveals its potential role in memory and learning (Zeisel et al. 2015). The Tac1 gene plays a role in regulating emotional responses, particularly in the mechanisms of depression (Quinn et al. 2013). Research on the Tcf7l2 gene further confirms its key role in thalamic neuronal functions, affecting neuronal firing patterns and behavioral choices (Lipiec et al. 2020). Knockout of the Zbtb20 gene impacts learning and memory, highlighting its role in hippocampal-dependent cognitive functions (Nielsen et al. 2014; Nagao et al. 2016). These findings emphasize how high-precision gene expression analysis techniques across different brain regions can deepen our understanding of how genes regulate brain development and disease processes on structural and functional levels.

To assess the impact of the microenvironmental deconstruction module on spatial domain identification accuracy, we conducted comparative experiments by replacing this module with other well-established models, including DestVI (Lopez et al. 2022), SPACEL*(Spoint), Cell2location (Kleshchevnikov et al. 2022), and Tangram (Biancalani et al. 2021). These models are widely recognized and applied in the field of deconvolution, providing robust benchmarks for our evaluation. Using the same data set and experimental conditions, we substituted the microenvironmental deconstruction module in the STMSC framework with each of the aforementioned models. The ARI was employed as the primary evaluation metric, measuring the consistency between clustering results and true labels in the spatial domain identification task. As shown in Figure 3D (purple bar graph on the right), the original microenvironmental deconstruction module in STMSC achieved the highest ARI score, outperforming the alternative models. This demonstrates the superior efficacy and applicability of our module in capturing spatial domain features and enhancing identification accuracy. In addition, we benchmarked performance against other deconvolution methods in terms of estimating cell type proportions using a simulation approach and found that STMSC exhibited comparable performance to existing methods for this specific task, albeit Cell2Location displayed the best performance as measured by a Pearson's correlation between predicted and observed cell type proportions (Cell2Location R = 0.867, STMSC R = 0.5125) (Supplemental Tables 2–6). These results highlight that different methods may be better suited to different analytical goals; for example, if precise cell type deconvolution is the primary objective, purpose-built methods such as Cell2Location may be preferable. Importantly, as detailed in the Methods section “Flexibility of the framework,” STMSC is designed as a modular and adaptable framework that allows users to incorporate the most appropriate deconvolution methods based on the specific characteristics and analytical goals of their data sets.

STMSC provides deeper insights into complex cancer environments through accurate spatial domain identification

Tumor heterogeneity presents substantial challenges in designing effective therapeutic strategies and is associated with poorer patient survival outcomes (Alizadeh et al. 2015; Nguyen et al. 2021). ST, by preserving spatial information crucial for reconstructing spatial distribution domains and distance-dependent interactions, offers unique advantages in anatomically characterizing and delineating intra-tumor heterogeneity. We applied STMSC to human breast cancer data acquired via the 10x Genomics Visium platform. Noted for its marked intra- and intertumoral heterogeneity, the data set was annotated into 20 spatial domains based on pathological features in the SEDR study (Xu et al. 2024), categorized into four morphotypes: ductal carcinoma in situ/lobular carcinoma in situ (DCIS/LCIS), healthy tissue, invasive ductal carcinoma (IDC), and tumor margins (Fig. 4A). First, in the task of spatial domain identification, using SEDR's annotations as a reference, Figure 4A illustrates the visualization of spatial domain recognition, demonstrating accurate and smooth identification for most regions. Figure 4H employs the ARI to quantify STMSC's performance in spatial domain identification, showing superior performance across most regions compared with alternative methodologies. Figure 4B illustrates the correspondence between different STMSC spatial regions and the manually annotated regions.

Figure 4.

Human breast cancer data set analysis results. (A) 2D visualization of STMSC spatial domain identification compared with true labels. (B) Alluvial plot showing correspondence between STMSC domains and manual annotations (ground truth). (C) Differential gene expression analysis for spatial domains 9 and 10. (D) Spatial expression visualization of CXCL14 and COX6C in spatial domains 9 and 10. (E) GO enrichment analysis with spatial domains 9 and 10 cross-referenced to each other. (F) GO enrichment analysis for spatial domains 0, 4, and 17. (G) GO enrichment analysis with spatial domains 2 and 6 cross-referenced to each other. (H) Comparison of ARI clustering evaluation metric for STMSC.

Open in new tabDownload PowerPointLink to figure

Second, we conducted differential gene expression (DEG) (Fig. 4C) and Gene Ontology (GO) enrichment analysis (Fig. 4E) on the identified spatial domains to validate the SEDR annotations and explore the heterogeneity between two spatial domains within the IDC region. First, we performed a differential gene expression analysis between spatial domains 9 and 10. The results revealed that CXCL14 was overexpressed in domain 9, promoting tumor cell invasion and migration, particularly in breast cancer, in which it enhances tumor growth and metastasis by secreting secondary cytokines. This finding aligns with the invasive characteristics of IDC (Hara and Tanegashima 2012). Elevated expression of CCND1 in domain 9 is recognized as a critical driver of IDC progression, as its overexpression accelerates the cell cycle, contributing to uncontrolled cell proliferation, a hallmark of IDC (Mercapide et al. 2002; Tajiri et al. 2014). In contrast, COX6C was highly expressed in domain 10, potentially facilitating tumor cells’ evasion of apoptosis through mitochondrial regulation, making it a promising biomarker or therapeutic target (Pareja et al. 2020; de Oliveira et al. 2023). The high expression of SLC39A6 in domain 10 has been linked to the progression of estrogen receptor-positive (ER⁺) IDC, suggesting its role as a potential marker for hormone therapy response (Thakkar et al. 2010). These findings underscore the consistency between the SEDR annotations and the biological characteristics of IDC across domains 9 and 10. We further visualized the spatial distribution of the aforementioned differentially expressed genes, CXCL14 and COX6C, as examples (Fig. 4D). The results revealed clear spatial patterns for both genes.

Furthermore, recognizing the complexity of the tumor microenvironment, even within the IDC region, heterogeneity can still arise. To probe the differences between domains 9 and 10, we performed a GO enrichment analysis using each domain as a reference for the other (Fig. 4E). The results demonstrated that domain 9 was enriched for GO terms related to extracellular matrix (ECM) organization and remodeling, processes essential for tumor cells to breach the basement membrane and invade adjacent tissues. Additional GO terms were associated with cell migration and motility, highlighting the enhanced invasive and metastatic potential of cells in this domain. In contrast, domain 10 was enriched for GO terms related to T cell activation and antigen receptor signaling, suggesting heightened immune activity in this area, potentially fostering a microenvironment conducive to immune cell infiltration and activation, which could impact tumor growth and progression. In summary, although both domains 9 and 10 belong to the IDC region, they exhibit heterogeneity. Domain 9 is characterized by increased cellular migration, invasion, and ECM remodeling, indicative of more aggressive and metastatic behavior, consistent with the invasive nature of IDC. In contrast, domain 10 displays enhanced immune-related activity, particularly in T cell activation and antigen receptor signaling, implying a more active immune response in the tumor microenvironment, potentially inhibiting invasion and metastasis. Therapeutically, domain 9 may benefit from strategies targeting invasion and metastasis, whereas domain 10 might be more responsive to immune-based therapies. Overall, this in-depth analysis of tumor-associated spatial domains not only validates the accuracy of the SEDR annotations but also uncovers the intraregional heterogeneity within IDC, providing novel insights into the complex biology of breast cancer and laying a robust foundation for future precision medicine research.

To further investigate the differences between distinct stages of breast cancer progression, we conducted GO enrichment analysis on spatial domains 2 and 6 (Fig. 4G). Using spatial domain 6 as a reference, GO analysis of spatial domain 2 revealed an enrichment of immune-related processes, particularly those associated with active immune responses and antigen presentation. This heightened immune activity in spatial domain 2 likely reflects the immune system's effort to detect and eliminate early-stage abnormal cells. Such immune engagement is characteristic of early breast cancer stages, such as DCIS and LCIS, in which the immune system plays a crucial role in responding to emerging tumor cells. Conversely, when spatial domain 2 was used as the reference, GO enrichment analysis of spatial domain 6 revealed an overrepresentation of pathways related to cell cycle dysregulation, uncontrolled cellular proliferation, and enhanced cell adhesion and migration. These processes suggest that cancer cells in spatial domain 6 have acquired invasive properties, consistent with the biological characteristics of IDC, in which cancer cells penetrate the basement membrane and invade surrounding tissues. These GO enrichment findings align with the known biological functions of spatial domains 2 and 6, further validating the spatial domain delineations. Moreover, they provide crucial molecular insights into the transition from in situ carcinomas, such as DCIS and LCIS, to invasive forms like IDC. This analysis underscores the potential molecular mechanisms driving breast cancer progression from localized, early-stage lesions to more aggressive, invasive disease.

Finally, it is important to note the complexity of tumor tissues, which, unlike brain cortex with clear and definite morphological boundaries, are highly heterogeneous and contain complex microenvironments. SEDR suggests that manual morphological annotations based solely on tumor pathology are insufficient to characterize this complexity. STMSC's recognition results for the upper right quadrant of pathological slices showed discrepancies; areas originally marked as healthy were identified as three distinct regions (spatial domains 0, 4, and 17). Detailed GO enrichment analysis (Fig. 4F) revealed that domain 0 is primarily involved in tissue remodeling and repair and domain 4 in immune response and signaling, and domain 17 reflects a combination of immune response, cell migration, and tissue repair. Given the relative spatial positions of domains 0, 4, and 17, these areas form what was artificially marked as a healthy region, with domain 17 on the outermost layer, domain 4 in the middle, and domain 0 at the center. We hypothesize that the immune activity in the outermost layer (domain 17) may be a response to cancer cell invasion, a notion further supported by consistent findings across multiple spatial domain identification methods (see Supplemental Section 1; Supplemental Fig. 1). This suggests that breast cancer may be attempting to invade the healthy regions composed of domains 0, 4, and 17, although these areas employ various mechanisms to protect themselves from harm. This insight underscores the efficacy of our data analysis and provides crucial guidance for future biological and clinical research.

Case study: application of STMSC in HER2-positive breast cancer data analysis

In our previous analysis of a single breast cancer slice, we demonstrated that STMSC enables accurate identification of spatial domains, offering deeper insights into the complexity of the tumor microenvironment. To further showcase the potential of STMSC in analyzing 3D tumor biology, we applied it to a HER2-positive breast cancer data set from patient E, which includes three adjacent tissue sections. Compared to single-slice analysis, multislice integration is critical for reconstructing the 3D spatial structure of tissues and obtaining a more comprehensive view of the tumor architecture. This approach enables a more continuous and spatially coherent interpretation of the tumor microenvironment.

By performing cross-slice integration using STMSC (Fig. 5B; Supplemental Fig. 2), we leveraged interslice continuity to enhance the robustness of spatial domain identification and recover biologically meaningful patterns along the depth axis of the tissue. Because only the first slice contains manual annotation (Fig. 5A), we hypothesized based on this slice's annotations that spatial domain 0 is likely immune infiltrate tissue, spatial domains 1 and 2 are likely invasive cancer tissue, and spatial domain 3 likely represents connective tissue. To further validate this hypothesis, we analyzed the enrichment or depletion of different cell types across space domains using the information obtained from STMSC (Fig. 5C) and conducted GO enrichment analysis (Fig. 5D). These analyses confirmed the biological relevance of STMSC-identified domains and provided deeper insights into the tumor microenvironment.

Figure 5.

HER2-positive breast cancer data set analysis results. (A) Pathologist's annotation for slice 1 of patient E. (B) Visualization of STMSC's spatial domain detection results. (C) Enrichment or depletion of different cell types in all regions. (D) Notable GO terms for regions.

Open in new tabDownload PowerPointLink to figure

First, through the analysis of enrichment or depletion of different cell types across all regions (Fig. 5D), we found enrichment of immune cells in spatial domain 0, including T cells CD4⁺ (Speiser et al. 2023), T cells CD8⁺ (Strioga et al. 2011), and NKT cells (Godfrey et al. 2004). The enrichment of plasmablasts further suggests that this region is a hotspot for immune cell infiltration, indicating a strong immune response (Mao et al. 2017; Wouters and Nelson 2018). In contrast, the depletion of cancer HER2 SC cells indicates that this region does not support tumor growth or expansion, possibly owing to immune-mediated inhibition of tumor cell proliferation (Lu et al. 2023). Thus, spatial domain 0 is identified as an immune-infiltrated region. Space domains 1 and 2, showing enrichment of cancer HER2 SC cells, suggest that these regions are closely associated with tumor growth, invasion, and metastasis, reflecting the characteristics of invasive cancer (Wang et al. 2024). The presence of mature luminal cells in these regions indicates epithelial development, differentiation, and remodeling, processes commonly associated with aggressive tumor progression (Shehata et al. 2012). Additionally, the enrichment of CXCL12-expressing endothelial cells suggests that these regions promote angiogenesis and microenvironmental changes, supporting tumor expansion and dissemination (Portella et al. 2021; Cambier et al. 2023). In contrast, spatial domain 3 shows enrichment of immune cells (e.g., macrophages, monocytes, and dendritic cells) as well as MSC iCAF-like and myCAF-like cancer-associated fibroblasts (CAFs). The coenrichment of immune cells and CAFs suggests a prominent presence of connective tissue cells, particularly fibroblasts (CAFs), which play a key role in tumor microenvironment remodeling and immune evasion. As the predominant cell type in connective tissue, CAFs are involved in tissue repair and reorganization of the tumor microenvironment (Maia et al. 2023; Tsoumakidou 2023). Therefore, spatial domain 3 is likely a connective tissue region that supports tumor progression and metastasis. Taken together, our analysis suggests that spatial domain 0 corresponds to an immune-infiltrated region, spatial domains 1 and 2 are consistent with invasive cancer tissue, and spatial domain 3 represents connective tissue.

To further validate these findings, we conducted GO enrichment analysis across the four spatial domains (Fig. 5D). For space domain 0, enriched pathways such as “B cell receptor signaling pathway” and “antigen receptor-mediated signaling pathway” are closely linked to the immune microenvironment in breast cancer, particularly in HER2-positive cases, in which immune activation plays a critical role in anti-HER2 therapy (Duro-Sánchez et al. 2023). The enrichment of the “B cell receptor signaling pathway” suggests active B cell function and antigen presentation in this region. In spatial domain 1, GO analysis revealed enrichment in “peptide antigen processing and MHC protein complex” pathways, indicating high immune activity in this region (Hewitt 2003). The “positive regulation of macrophage migration” further suggests that this domain is a major site for immune cell infiltration, particularly macrophages, and is associated with immune escape and tumor progression (Huang et al. 2024). GO analysis of spatial domain 2 showed enrichment in pathways related to cell migration and apoptosis regulation, both of which are closely associated with the metastatic potential of HER2-positive breast cancer. The notable enrichment of apoptosis-related processes suggests that this region may contribute to tumor cell survival and resistance to cell death (Wani et al. 2023; He et al. 2024). For space domain 3, GO enrichment analysis again highlighted immune-related pathways such as “B cell receptor signaling pathway” and “antigen receptor-mediated signaling pathway,” emphasizing its involvement in immune response. In HER2-positive breast cancer, immune evasion is a critical factor in tumor progression (Gil Del Alcazar et al. 2020). The enrichment of “immunoglobulin-mediated immune responses” further supports the idea that this region may contribute to immune evasion or responses to immunotherapy (Wang et al. 2020a; Zhang and Kim 2024). Overall, the GO enrichment analysis of the four space domains is consistent with the conclusions drawn from the cell type enrichment/depletion analysis, further validating the functional features of these domains and providing valuable biological insights into the tumor microenvironment.

Overall, this case study demonstrates the potential of STMSC to reconstruct 3D tissue architecture through multislice integration, enabling a comprehensive and biologically meaningful interpretation of complex tumor structures.

Preprocessing

To ensure reproducibility in data analysis, STMSC implemented routine preprocessing of both ST and SC data sets using the SCANPY package (Wolf et al. 2018). The preprocessing protocol involves several key steps: Initially, mitochondrial-encoded genes are removed, and gene identifiers are standardized. Common genes shared by the ST and SC data sets are retained and normalized to mitigate systemic sample biases. Subsequently, STMSC selects the top 500 highly variable genes (HVGs) from the SC data set, which exhibit expression variability and provide crucial biological insights. Based on these HVGs, the average expression values for each cell type are calculated, forming the final data set, X_sc. In this data set, rows represent specific cell types, columns represent genes, and $x_{ij}\in {\mathbf{X}}_{sc}$ denotes the average expression level of a gene j within a particular cell type i.

In extracting information from H&E images, STMSC employed a pixel coordinate-based mapping strategy to align specific spots with H&E images, in which the average RGB value (r_m, g_m, b_m) of the pixels in the 50 × 50 windows around the coordinate point is used to determine specific pixel coordinates.

Slice registration

Traditional spatial coordinate-based alignment algorithms, which typically rely on defining edge spots for slice alignment, often prove inadequate when dealing with samples characterized by irregular tissue margins or complex internal structures (Wang et al. 2023). Additionally, the reliance on point-to-point matching algorithms based solely on coordinate data renders these methods particularly vulnerable to interference from outliers, errors that frequently stem from the inherent limitations of sequencing technologies. Such inaccuracies can compromise alignment precision. To overcome these challenges, the STMSC model introduces a multislice alignment module that integrates spatial coordinates (x_pc, y_pc) with histological features (r_m, g_m, b_m) from H&E images. This approach initially utilizes spatial coordinate data from multiple slices as input and applies the iterative closest point (ICP) (Arun et al. 1987) algorithm for initial alignment. The method then enhances this alignment by incorporating histological image features, treating them as a third dimension alongside the coordinate data (Hu et al. 2021), and refines the alignment using the ICP algorithm again. In contrast to atlas-based alignment frameworks, such as the Allen Common Coordinate Framework (CCFv3) (Wang et al. 2020b), a high-resolution 3D brain atlas commonly used for standardized registration in brain studies, STMSC adopts a template-free alignment strategy. This design choice reflects the practical need to accommodate human and cancer tissue samples, which are incompatible with CCFv3 owing to species- and organ-specific differences. By leveraging intrinsic spatial and histological features, STMSC offers a generalizable framework applicable across a wide range of tissue types, especially in contexts in which standardized anatomical references are unavailable.

To incorporate the RGB image information from H&E staining as 3D information, we convert the 3D RGB data into a single-dimensional value through a weighted sum. The conversion formula is as follows:

The multislice alignment step takes into account more complex biological characteristics and histological differences, thereby identifying and correcting regions that cannot be accurately aligned based solely on coordinate information. This ensures that matched points are not only spatially proximate but also histologically similar. The final alignment result is

Microenvironment deconstruction

By integrating SC data, STMSC effectively deconstructs the spot-level microenvironments, transforming coarse-grained information into a fine-grained representation at the cellular level. This enhanced granularity is achieved through a nonaugmented contrastive learning model that accurately deconvolutes the distribution of cell types within each spot. As a result, STMSC effectively reconstructs the cellular landscape, providing a more detailed and accurate view of the tissue microenvironment. Specifically, STMSC constructed a mapping matrix M of dimensions N_celltype × N_spot to map the probability of cell type within a spot. Each element m_ij in this matrix signifies the probability of mapping the ith cell type to spot j. The mapping matrix M is learned through the equation

The contrastive loss component emphasizes maximizing the similarity of positive relations between spatially adjacent spots while minimizing the similarity of negative relations between nonadjacent spots. The reconstruction loss component ensures alignment between predicted and reconstructed gene expression profiles. Regulatory factors α and β (set to α, β ∈ [0, 1], α + β = 1) are introduced to finely balance the weights of contrastive and reconstruction losses, optimizing overall model performance.

Construction and correction of 3D neighborhood graphs

The STMSC model leverages the aligned spatial positional information by employing Euclidean distance to compute the distances between points, thereby generating a distance matrix D. The values within D represent the strength of relationships between two spots. A suitable threshold s is set for different data sets, such that if d_ij < s, the spots i and j are considered physically close.

The microenvironment deconstruction module in STMSC constructs the cellular landscape, which then serves as input for the precorrection mechanism. These fine-grained data are crucial for correcting graph structures, leading to a more detailed understanding of the relationships between spots. Enhanced granularity improves the identification of spatial domains, ultimately advancing our comprehension of tissue architecture and function. To achieve this, the model incorporates fine-grained information M obtained from microenvironment deconstruction for correction. This allows for the consideration of cellular composition and heterogeneous states within microenvironments. For example, if physically proximate spots exhibit differences in cellular states or composition, the connections in D may be downweighed or redefined to accurately reflect actual biological differences.

Based on the corrected distance matrix D^*, an initial 3D neighborhood graph $\mathbf{A}\left(N_{spot}\times N_{spot}\right)$ is constructed, in which nodes represent individual spots and edge weights are determined by the threshold s. If $d_{ij}^{\ast }<s$, it indicates that spot i and spot j are neighbors in 3D space, and hence, a_ij = 1; otherwise, a_ij = 0.

Network structures

Within the STMSC framework, an autoencoder is engineered to learn complex biological information derived from 3D biological landscapes. The encoding network within the autoencoder is designated as F_en(X_st, A). To enhance the model's capabilities, we incorporate a graph attention network (Veličković et al. 2018) that consists of graph attention layers and dense layers. The attention mechanism enables dynamic weighting of the influence of one spot on another, based on the spatial and molecular context of each spot. This integration enhances the model's capacity to interpret and capture complex spatial relationships among spots. Compared to CGNs, GAT enables adaptive weighting of neighboring nodes, which is better suited for the heterogeneous spatial context of transcriptomic data. In addition, it offers improved scalability and generalization over transformer-based models, especially on small or noisy data sets. The input to the encoding network F_en includes a 3D global adjacency matrix, A, and preprocessed ST gene expressions, X_st. The output dimension of the graph attention layer is set at 512. Subsequently, the dense layer in F_en receives these 512-dimensional hidden vectors and outputs 128-dimensional latent representations, Z. Similarly, the decoder F_de(Z) of the autoencoder in the STMSC framework utilizes a construction analogous to the encoder, featuring a single-layer dense network and a graph attention network with a hidden dimensionality of 512. This network takes the latent representations, Z, as input and ultimately reconstructs the gene expression matrix through a dense layer.

The loss function of this GAE is denoted as

Clusters

The output Z from the hidden layer of the GAE serves as input for clustering algorithms tailored to delineate biologically meaningful spatial domains. The choice of clustering strategy depends on the availability of prior information regarding the number of spatial domains. Specifically, for data sets with known domain counts, a Gaussian mixture model (GMM) (McNicholas and Murphy 2010) is employed for clustering. Conversely, in the absence of such information, the Louvain algorithm (Blondel et al. 2008) is utilized. The accuracy of clustering outcomes is assessed using the ARI (Hubert and Arabie 1985):

Hyperparameter settings

To ensure a fair and rigorous comparison, the hyperparameter settings for each benchmarked method were carefully selected based on the guidelines and strategies provided in the official tutorials and code repositories by the respective method authors. By adhering closely to these recommended settings, we aimed to minimize bias and maximize the validity of our experimental results, thereby ensuring that each method was evaluated under optimal and comparable conditions:

Flexibility of the framework

STMSC is designed as a flexible framework that can be adapted to different data characteristics and analysis needs. Rather than being constrained to a specific deconvolution module, STMSC allows users to integrate suitable deconvolution methods based on the properties of their data sets. This adaptability ensures that the framework remains applicable across diverse ST studies, accommodating variations in tissue complexity, and technological platforms. By enabling modular customization, STMSC encourages further exploration and refinement, making it a versatile tool for spatial domain identification and downstream biological analysis.

Data sets

The ST and scRNA-seq data sets used in this study are publicly accessible. The DLPFC data set can be obtained from http://research.libd.org/spatialLIBD/ (Maynard et al. 2021). The mouse coronal brain data set from the 10x Visium platform is available at https://squidpy.readthedocs.io/en/latest/notebooks/tutorials/tutorial_visium_hne.html (Palla et al. 2022), whereas the human breast cancer data set from the 10x Visium platform can be accessed at https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0. For the scRNA-seq data, the human DLPFC data set, analyzed using the 10x Genomics Chromium platform, is available from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE144136 (Nagy et al. 2020). The mouse brain data set, also analyzed using the 10x Genomics Chromium platform, can be found in the ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress) under accession number E-MTAB-11115 (Kleshchevnikov et al. 2022). Additionally, the SC transcriptomic data for human breast cancer is accessible at https://singlecell.broadinstitute.org/single_cell/study/SCP1039 (Wu et al. 2021). The HER2-positive breast tumor data set (Andersson et al. 2021) was profiled by ST platform (https://doi.org/10.5281/zenodo.4751624); the HER2-positive breast tumor data set (Wu et al. 2021) was profiled by 10x Genomics Chromium platform37 (GEO; GSE176078).

Software availability

The STMSC framework is available as an open-source project on GitHub (https://github.com/bliulab/STMSC) and as Supplemental Code. It is also freely accessible via the project website (http://bliulab.net/STMSC).