Spatially Resolved Transcriptomics Data Clustering with Tailored Spatial-scale Modulation

1 Introduction

Spatial transcriptomics technology [Song and Su, 2021; Zhai et al., 2023; Rao et al., 2021] provides a novel approach to studying the functions and interactions of cells within tissues. Offering quantitative gene expression data alongside the spatial location of cells overcomes the limitations of single-cell sequencing [Mrabah et al., 2023], particularly the loss of spatial information. In recent years, spatial transcriptomics has been widely applied across various fields, including embryonic development [Rossant and Tam, 2022], disease research [Moore et al., 2024], and plant science [van Dijk et al., 2021].

∗Corresponding author.

Accurate identification of spatial domains within transcriptomic data enables a deeper understanding of the distribution patterns of cell populations in tissues and their biological functions at specific spatial locations. Early methods for spatial domain identification primarily relied on gene expression data for clustering analysis, such as K-means [Ahmed et al., 2020], Louvain [Zhang et al., 2021], and Seurat [Pereira et al., 2021] algorithms. However, these methods neglected spatial relationships between cells, often leading to discontinuous domain identification, where adjacent cells were mistakenly assigned to separate domains.

Consequently, researchers have increasingly explored deep clustering methods integrating gene expression data with spatial location information. For example, stGCL [He et al.,

2024] constructs a heterogeneous graph coordinating module that maintains view independence during the fusion process. stMMR [Zhang et al., 2024] introduces a multimodal geometric method to effectively fuse histological, gene expression, and spatial location information. SCGDL [Liu et al., 2023a] applies a Bayesian Gaussian mixture model to generate spatial domains. Recently, many methods based on Graph Neural Networks (GNNs) have also been proposed, such as MAFN [Zhu et al., 2024], SpaGCN [Hu et al., 2021], and Spatial-MGCN [Wang et al., 2023]. However, these methods construct spatial graphs using a conventional approach with a static radius based on spatial coordinates, as shown in Fig. 1 A. For the central node, considering all nodes within the circle as neighbors, including those from other clusters, hinders the accurate identification of spatial domains and complicates the precise partitioning of boundary nodes in clusters. The radius used to construct the graph is crucial: if it is too small, fails to capture adequate spatial neighbors; and too large, it introduces noise into the spatial domain.

To address this issue, we propose a novel spatially resolved transcriptomics data clustering network (TSstc) with tailored spatial-scale modulation, as illustrated in Fig. 2. Specifically, we employ a tailored spatial-scale modulation approach that incrementally constructs different spatial graphs Avspa as the spatial domain radius expands. Additionally, we introduce a Spatiality-Aware Sampling (SAS) strategy to aggregate node representations by accounting for the spatial dependencies between spots. Next, we use graph convolutional network (GCN) encoders to learn gene embeddings Zgene from gene graph and multiple spatial embeddings Zvspa from spatial graphs. The Cross-view Fusion module is then applied to integrate the gene embedding and spatial embedding features for Zf . During training, we incorporate cross-view correlation-based tailored spatial regularization constraints across spatial embeddings Zvspa at different scales to remove noisy neighbors while preserving high-quality neighbor relationships. Finally, we utilize a zero-inflated negative binomial (ZINB) model to reconstruct the original expression matrix, capturing the global probability distribution of gene expression profiles. Extensive experimental results demonstrate that our approach outperforms existing state-of-the-art methods in various spatial transcriptomic tasks and related downstream applications.

The main contributions of this work are summarized as follows:

• We propose a novel network for clustering spatially resolved transcriptomics data, termed TSstc. To our knowledge, we are the first to explore spatial graph construction scales and address the challenge of identifying boundary nodes within clusters.
• We design a tailored spatial-scale modulation approach that incrementally constructs spatial graphs and employs the Spatiality-Aware Sampling (SAS) strategy to aggregate node representations. Additionally, we implement cross-view correlation-based tailored spatial regularization constraints to preserve high-quality neighbor relationships across spatial embeddings at different scales.
• We conducted experiments on three benchmark datasets

to validate the effectiveness of our proposed method. The results demonstrate the superiority of TSstc compared to state-of-the-art methods.

2 RelatedWork

2.1 Spatial Transcriptomics Techniques

Spatial transcriptomics (ST) techniques [Khan et al., 2024; Yang et al., 2024c], widely used in fields such as embryonic development [Rossant and Tam, 2022], disease research [Moore et al., 2024], and plant science [van Dijk et al., 2021], provide insights that traditional transcriptomics cannot offer due to its lack of spatial resolution. ST enables the capture of quantitative gene expression data while preserving the spatial locations of cells, thereby overcoming the spatial information loss inherent in single-cell sequencing. As a result, significant efforts have been made to develop deep learning networks that utilize spatial domains. Hu et al. [Hu et al., 2021] proposed a graph convolution-based network that aggregates gene expression of each spot and its neighboring spots to identify coherent spatial domains. Wang et al.[Wang et al., 2023] introduced an attention mechanism to combine gene expression and spatial embeddings, followed by spatial regularization to constrain the integrated features. However, these methods typically map aggregated features into a unified representation, overlooking gene heterogeneity. Xiao et al. [He et al., 2024] constructed a heterogeneous graph coordinating module that maintains view independence during the fusion process. Zhu et al.[Zhu et al., 2024] designed a Cross-view Fusion Module to dynamically integrate multiple features and employed a CCR strategy to reduce irrelevant information. Recently, Graph Neural Networks (GNNs) have shown remarkable effectiveness in handling non-Euclidean structural data, leading to their growing application in spatial transcriptomic clustering. Notable examples include SEDR [Xu et al., 2024], DeepST [Xu et al., 2022], CCST [Li et al., 2022], GraphST [Long et al., 2023], and others. These approaches leverage graph deep learning to accurately model the spatial distribution and interactions of gene expression within cells, thereby improving the precision of spatial domain delineation.

2.2 Deep Learning-basedMulti-viewClustering

In recent years, deep learning-based multi-view clustering methods have attracted significant attention. These methods can be broadly categorized into three groups: generative adversarial network-based (GAN) [Yang et al., 2024b; Yang et al., 2022], variational autoencoder-based (VAE) [Li et al., 2023; Yuan et al., 2023], and graph convolutional network-based (GCN) [Xiao et al., 2025]. Specifically, GANs generate realistic view data samples to improve clustering performance. Li et al. [Li et al., 2019] used an adversarial module to capture the multi-view data distribution and disentangle the latent space. Meanwhile, VAEs map input data into a probabilistic distribution in the latent space, from which latent variables are sampled. Xu et al. [Xu et al., 2021] have defined view-common and view-specific variables within a generative model, extracting a common clustering factor through an approximation of the discrete Gumbel- Softmax distribution. Additionally, recent literature also benefits from the strong representation power of graph convolutional networks, showing a strong tendency to utilize GCNs to further improve clustering performance. Zhou et al. [Zhou et al., 2024] developed a self-encoding framework that incorporates view-specific GCN networks to learn high-level representations of each view's data features while considering prognosis prediction results. In our proposed method, TSstc, we aim to explore how to integrate gene information and spatial relationships to improve clustering performance.

3 Methodology

In this section, we first provide an overview of the TSstc model. We then describe the individual modules of the proposed model in detail, followed by an explanation of the training loss function.

3.1 Overview of TSstc

Our proposed TSstc model consists of four main components: the Graph Construction module, the Spatiality-Aware Sampling (SAS) module, the Feature Extraction module, and the Spatial Domain Constraints module.

The TSstc model is trained on spatial transcriptomic datasets, which include the gene feature matrix Xgene ∈ RN×D and the 2D spatial coordinates Sc ∈ RN×2 , where N denotes the number of cells and D denotes the dimension of each cell. We construct the gene graph Agene ∈ RN×N and different spatial graphs {Avspa }nv=1 ∈ RN×N incrementally as the radius of the spatial domain expands. Then, we aggregate node representations to obtain different spatial feature matrixes {Xvspa }nv=1 ∈ RN×D by the SAS module and we utilize GCN to learn the latent representations Zgene ∈ RN×d and Zvspa ∈ RN×d . To obtain more comprehensive crossview features Zf , TSstc uses an attention fusion module to integrate the gene and spatial features. Furthermore, we use the Spatial Domain Constraints module to preserve high-quality spatial neighbor relationships and capture the global probability distribution of gene expression profiles.

3.2 Graph ConstructionModule

Compared to directly constructing spatial graphs with a fixed radius based on spatial coordinates Sc, we employ a tailored spatial-scale modulation approach to construct multiple spatial graphs. For each dataset, we analyze the spatial distribution of spots and incrementally construct spatial graphs Avspa as the Euclidean distance-based radius Rv increases, thereby evaluating the connectivity between two specified spots as follows:

where Lvij represents the Euclidean distance of spot i and j, calculated based on the spatial coordinates Sc.

Then we measure the cosine similarity between two points i and j by calculating the similarity of gene expression vectors Xgenei and Xgenej . We then construct the K-NN [Guo et al., 2003] neighbor graph, denoted as Agene, using a binary classification system, where:

We utilize the tailored spatial-scale modulation approach to construct various spatial graphs Avspa that identify spatial domains of spots. Additionally, the gene graph Agene is constructed to capture the similarity between spots.

3.3 Spatiality-aware Sampling (SAS) Module

Understanding the spatial domain of the target spot and sampling its spatial neighborhood is crucial for spatially resolved transcriptomics data. We propose the Spatiality-Aware Sampling (SAS) module and combine the various spatial graphs Avspa to perform an aggregation operation on the features of the target spot and its spatial neighborhood, yielding the spatial feature representation Xvspa before training, where:

where Ni represents the spatial neighborhood of spot i, and the aggregation function can take various forms, such as mean, max, sum, attention, and ensemble, which learn the spatial feature representations. In this case, we choose the mean operation.

3.4 Feature ExtractionModule

In general, spatially resolved transcriptomics data are prone to noise and redundant information. To filter out this noise and retain more discriminative features in the latent space, we use the GCNgene encoder and GCNspa(v) encoders [Yang et al., 2024a; Qin et al., 2024] to learn specific features of the gene graph and spatial graph, respectively. The propagation rule is defined as follows:

We simplify the gene and spatial normalized adjacency matrix as ˜A, which is given by ˜A = (D)- 1 2 (I + A)(D)- 1 2 , where Dii = ∑j (Aij + Iij), and the identity matrix I ∈ RN×N ensures self-loops for each spot. The matrices W(l) 1 and Wv(l) 2 are the weights for the l-th layer and σ (·) is the Relu activation function.

In this subsection, we introduce the Cross-view Attentioninduced Feature Fusion module, which allocates variable weights to effectively integrate gene expression and spatial embedding features. Specifically, we integrate the latent representations Zgene ∈ RN×d and {Zvspa }nv=1 ∈ RN×d to obtain the global cross-view feature Zf ∈ RN×d .

To incorporate the (v + 1)-th view information, we concatenate the d-dimensional feature vectors of Zvspa and Zgene along the feature dimension to form Zc ∈ RN×(d×(v+1)) . An attention weight matrix H ∈ RN×(v+1) is then generated by the network as follows:

where ϕ(·) denotes the Softmax function, δ (·) is the Tanh activation function, andW1, W2, b1, and b2 are parameters of the fusion network.

Next, a broadcasting mechanism is applied to ensure dimensional compatibility, followed by the application of attention weights to each feature vector. Finally, the weighted feature vectors are summed to yield the fused feature vector Zf ∈ RN×d as follows:

where Hv denotes the v-th column of the attention weight matrix H and Hlast is its last column.

3.5 Spatial Domain Constraints Module

Raw spatial transcriptomics data include spots with spatial relationships and location information, making it crucial to preserve these relationships in the latent space. SpatialMGCN [Wang et al., 2023] and MAFN [Zhu et al., 2024] use spatial graph-based regularization loss to enforce embedding constraints. However, if the spatial graph contains noisy nodes or connections between nodes from different clusters, the regularization loss hinders the learning of highquality embedding features. To address this, we propose a cross-view correlation-based tailored spatial regularization constraint loss to preserve high-quality neighbor relationships across spatial embeddings at different scales. We begin by evaluating the sample correlation Sv between the gene features Zgene and each spatial feature Zvspa , as follows:

where Svij represents the cosine similarity between the i-th node embedding in the v-th spatial view and j-th node embedding in the gene view. For each spot in the spatial embedding feature Zvspa , we select its neighboring nodes with high expression values based on the correlation matrix Sv to refine the spatial graph. This process removes some connections that may originate from different clusters or contain noise, with low expression values. We then apply regularization loss to the spatial embedding features, as shown below:

where Rvi denotes the set of spatial neighbors for spot i in the v-th spatial embedding. It is generally accepted that data spots in close physical proximity should also be close in latent space, while non-adjacent spots should be distant.

As done in previous methods [He et al., 2024; Zhu et al., 2024], we use the ZINB decoder module [Wang et al., 2023] to reconstruct the gene expression matrix, assuming that the gene expression matrix Xgene and the spatial expression matrix Xspa follow a ZINB distribution. For simplicity, we denote them as X in the following equations:

where pnb (Xij | bi) is given by:

Here, πij and vij are the zero-inflation and mean parameters, respectively, while θij and bi represent the discretization of the decoder output and the bias vector. The latent representation Zgene and Zvspa serve as input to the ZINB decoder, yielding estimates for these parameters. The negative loglikelihood of the ZINB distribution is used as the reconstruction loss:

3.6 The Total Loss Function

Our total loss function is defined as follows:

The parameter α represents the weighting factor used to balance the influences of LR and LC.

4 Experiments

This section assesses the performance of our approach on three benchmark datasets and compares its clustering effectiveness with that of eight state-of-the-art methods for spatial transcriptomics.

4.1 Datasets

We evaluate our method on three widely used spatial transcriptomics datasets. The first dataset is the LIBD human dorsolateral prefrontal cortex DLPFC [Maynard et al., 2021], comprising 8 slices each containing 33,538 genes and 5-7 regions. The second dataset is the 10x Visium dataset of human breast cancer [Maynard et al., 2021] which contains 20 regions and 36,601 genes. The third dataset is mouse brain anterior tissue [Buache et al., 2011], annotated with 52 regions and containing 32,285 genes. Preprocessing for these datasets follows a method similar to that described in MAFN [Zhu et al., 2024]. We exclude spots outside the primary tissue areas and use the SCANPY [Wolf et al., 2018] toolkit to filter out genes with minimal expression or variance, selecting the top 3,000 genes with the highest variability [Zeng et al., 2023]. Finally, gene normalization is performed by applying a scaling factor.

4.2 ComparedMethods

To validate the effectiveness of our method, we compare it against eight approaches: SCANPY [Wolf et al., 2018], SpaGCN [Hu et al., 2021], DeepST [Xu et al., 2022], GraphST [Long et al., 2023], SCGDL [Liu et al., 2023b], Spatial-MGCN [Wang et al., 2023], stLearn [Pham et al., 2023], stCluster [Wang et al., 2024]. All comparative experiments are conducted using default parameters.

4.3 Training Details andMetrics

We conduct experiments to evaluate the proposed TSstc on the PyTorch platform using a single NVIDIA GeForce RTX 3090. Identifying the optimal spatial domain is crucial for different datasets. First, we analyze the spatial distribution of spots and incrementally construct spatial graphs as the Euclidean distance-based radius R increases. Next, we input these spatial graphs into the GCN and train TSstc on all benchmark datasets for at least 150 iterations until convergence. We apply a weight decay of 0.0005 and use the Adam optimizer with a learning rate of 1e-3 to optimize our model. The parameter α is set to 10. To validate the effectiveness of clustering, we use the Adjusted Rand Index (ARI) [Meil˘a, 2007] and Normalized Mutual Information (NMI) [Knops et al., 2006] as evaluation metrics.

4.4 Clustering Results and Analysis

Table 1 summarizes the performance of TSstc in terms of ARI compared to other methods across ten datasets, with NMI provided in the supplementary materials. From this table, we observe the following:

(1) In the DLPFC dataset, TSstc outperforms nearly all methods across 8 slices, achieving 4.32% higher ARI scores in 151507 slice compared to stCluster [Wang et al., 2024], which uses graph contrastive learning. Several methods show significant drops below 40% on specific slices, reducing their reliability in scientific and clinical settings. In contrast, TSstc outperforms the average of competing methods by 23%, demonstrating the effectiveness of the tailored spatial-scale modulation approach.

(2) In the Human Breast Cancer dataset, our tailored spatial-domain modulation improves ARI by 0.4% and NMI by 1.94% compared to the best-performing model, Spatial-MGCN [Wang et al., 2023]. We also detect the spatial domain, and TSstc clusters similar spots effectively, leading to well-defined boundaries, as shown in Fig.3. In contrast, DeepST[Xu et al., 2022] and GraphST [Long et al., 2023] fail to preserve spatial coherence, despite their strong local representation capabilities. By using the Spatiality-Aware Sampling (SAS) mechanism, TSstc integrates spatial information at multiple scales, ensuring a biologically meaningful connection between neighboring spots based on expression patterns. In the Mouse Brain Anterior dataset, Tab. 1 compares TSstc with other methods, and the supplementary materials show the spatial domain detection results. TSstc achieves the highest ARI of 44.82%, outperforming other methods, while stLearn [Pham et al., 2023], stCluster [Wang et al., 2024], and SpaGCN [Hu et al., 2021] report lower ARI scores of 38.42%, 36.88%, and 31.65%, respectively.

4.5 Ablation Studies

This section presents an ablation study to examine the impact of each component in our TSstc model. The ablation experiments are categorized into four sections, and Tab. 2 shows ARI performance across ten real datasets, with NMI results provided in the supplementary materials.

Effect of LR As shown in Tab. 2, w/o LR refers to the removal of the spatial domain constraints module. In the 151672 slice, this results in a 27% decrease in performance compared to TSstc, emphasizing the importance of spatial domain constraints in preserving spatial relationships and location details and integrating LR into the spatial embedding features ensures the retention of high-quality spatial neighborhood information.

Effect of LC w/o LC represents the removal of the ZINB decoder module. Analysis of the third row in Tab. 2 indicates that incorporating LC improves the average performance by over 9.9%, further validating the effectiveness of the ZINB decoder. This module reconstructs the gene expression matrix, mitigating noise issues caused by the discrete and sparse nature of spatial transcriptomics data.

Effect of SAS module w/o SAS refers to the omission of the Spatiality-Aware Sampling (SAS) module. In the Human Breast Cancer dataset, TSstc achieves approximately 7% higher ARI scores compared to w/o SAS, demonstrating the module's effectiveness. By leveraging the SAS mechanism, TSstc integrates spatial information across multiple scales, ensuring a biologically meaningful association between neighboring spots based on expression patterns.

Effect of TSM module To evaluate the effectiveness of our tailored spatial-scale modulation (TSM) module, we replace it with a static radius, following the approach used in Spatial-MGCN [Wang et al., 2023]. The results show that w/o TSM leads to a 4% lower ARI in the Mouse Brain Anterior dataset, with consistently lower ARI scores across other datasets. Furthermore, we evaluate the effect of different spatial graph combinations across all datasets, as shown in Fig. 4. For each dataset, spatial graphs are incrementally constructed as the Euclidean distance-based radius R increases, facilitating the identification of the optimal spatial graph combination. For the 151672 slice, the radius combination {2, 3, 4, 5} yields spatial graphs that achieve approximately 4% higher ARI scores and 3% higher NMI scores compared to the second-best combination. For the 151507 slice, the radius combination {2, 3, 4} achieved the best scores. Overall, these findings highlight the importance of the TSM module in identifying the optimal spatial domain, effectively capturing spatial relationships, and enhancing the understanding of cellular survival patterns.

4.6 Parameter Sensitivity Analysis

This section investigates the impact of the parameter α on clustering performance by varying its values within the set {0.01, 0.1, 1, 10, 100} across the DLPFC 8 slices, Human Breast Cancer, and Mouse Brain Anterior datasets. α serves as a weighting factor that balances the influences of LR and LC. As shown in Fig. 5, in the 151507 slice, the ARI increases gradually with α, while in the Mouse Brain Anterior datasets, it decreases or remains constant. The Human Breast Cancer dataset shows minimal performance changes, indicating that its clustering performance is less sensitive to variations in α. In our experiments, the parameter α is fixed at 10.

4.7 DownstreamApplications

We explore tumor tissue heterogeneity by examining the top 10 differentially expressed genes (DEGs) in the Healthy1 and DCIS/LCIS5 clusters, as shown in Fig. 6. The identified gene expression variations highlight regional molecular differences, underscoring TSstc's ability to capture the complex heterogeneity of Human Breast Cancer tissues. We also conduct several downstream applications, including UMAP visualization and gene imputation. Detailed results are provided in the supplementary materials.

5 Conclusion

In this paper, we propose a novel spatially resolved transcriptomics data clustering network (TSstc) to accurately identify spatial domains and precisely partition the boundary nodes within clusters. Specifically, we introduce a tailored spatial-scale modulation approach that incrementally constructs spatial graphs as the spatial domain radius expands, and a Spatiality-Aware Sampling (SAS) strategy is proposed to aggregate node representations. We incorporate cross-view correlation-based tailored spatial regularization constraints and employ a zero-inflated negative binomial model to capture the global probability distribution. Extensive experiments demonstrate that our method outperforms state-of-the-art clustering approaches and benefits related to downstream applications. For future work, we plan to explore the spatial distribution patterns of spots to refine clustering performance.

Acknowledgments

The work was supported in part by the National Natural Science Foundation of China under grants 62476258 and 62325604, and in part by the Natural Science Foundation of Hubei Province under grant 2025AFA113, and in part by the Fundamental Research Funds for National Universities, China University of Geosciences (Wuhan) under grant 2024XLB6.