AE-GCN is an integrative scheme that incorporates the AE and GCN learning processes, enabling tasks of spatial domain detection and data denoising (Figure 1).

Given the original expression X 0 ∈ R M × N (where M and N , respectively, denote the number of genes and spots) and spatial coordinates, the spatially neighboring network A ∈ R N × N and the enhanced expression data X ∈ R M × N are computed as the input of the combined learning process (see Methods). On the enhanced expression data X , AE-GCN employs AE to learn the low-dimensional representation (i.e., AE-specific representation H 1 l , l = 1 , … , L , where L is the number of total layers in AE) in each layer. With the spatially neighboring network A , AE-GCN utilizes GCN to learn the graph-constrained representation in each layer (i.e., GCN-specific representation H 2 l , l = 1 , … , B , where B is the number of layers in GCN or the encoder of AE). Then, AE-GCN transfers the AE-specific representations from the encoder to the corresponding GCN layer, thus generating the combined representation Y . Additionally, AE-GCN proposes a clustering-aware contrastive module to make the combined representation appropriate for spatial clustering.

When the learning process reaches convergence, the low-dimensional representation (i.e., Y ) of the last layer and the reconstructed expression data (i.e., X ′ ) can be used for downstream analytical tasks. The optimal representation enables AE-GCN to identify spatial domains interoperating with the Leiden method (Traag et al., 2019). The reconstructed expression data serve as the denoised profile, which overcomes the sparsity of the original data to improve differentially expressed gene identification (see Methods).

We employ AE to learn the useful representations from the expression data itself and assume that there are B layers in the encoder and L − B layers in the decoder. Specifically, the learned l th layer representation, H 1 l , can be obtained as follows: H 1 l = ϕ l W l H 1 l − 1 + b l , l = 1 , … , B (4) H 1 l = ψ l W l H 1 l − 1 + b l , l = B + 1 , … , L (5) where ϕ l and ψ l are the activation functions of the l th layer in the respective encoder and decoder. W l and b l are the weight matrix and reconstruction error in the l th layer, respectively. For convenience, we denote the enhanced expression data X as H 1 0 .

The output (i.e., X ′ = H 1 L ) of the decoder part is obtained through the reconstruction of the input data (i.e., X ) by minimizing the following loss function: L r e s = X − X ′ F 2 (6)

AE-specific representations, e.g., H 1 1 , H 1 2 , ⋯ , H 1 L , can denoise data itself and extract valuable information from the data itself, which can effectively reflect expression variation but cannot guarantee the spatial smoothness of the identified domains. GCN can model the spatial structural dependency between spots, which is beneficial to improving the spatial smoothness of the identified domains. Thus, we then transfer AE-specific representations in the encoder into GCN-specific representations and use the GCN module to propagate these AE-specific representations for capturing a more complete and powerful representation. Thus, the GCN-learnable representations can accommodate two different kinds of information: gene expression values and spatial neighborhood structure. The representation learned by the l th layer of GCN, H 2 l , can be obtained as follows: H 2 l = 1 − μ ϕ l W l H 2 l − 1 D ∼ − 1 2 A ∼ D ∼ − 1 2 + μ   H 1 l D ∼ − 1 2 A ∼ D ∼ − 1 2 (7) where I denotes the identity diagonal matrix. A ∼ = A + I and D ∼ i i = ∑ j A ∼ i j . D ∼ − 1 2 A ∼ D ∼ − 1 2 is the normalized adjacency matrix. μ is the balance coefficient and is often uniformly set to 0.5. Note that GCN and AE share weights.

Note that we denote the representation (i.e., H 2 B ) of the last GCN layer as Y . The input of the first layer GCN can be obtained from the enhanced expression data X : H 2 0 = ϕ l W 0 X D ∼ − 1 2 A ∼ D ∼ − 1 2 (8)

Although we have incorporated the encoder of AE into the neural network architecture of GCN to obtain the combined representation, this representation cannot be directly applied to the clustering problem. Herein, we propose a clustering-aware contrastive module to unify these two different deep learning models for effective spatial clustering.

Specifically, we use student’s t-distribution to measure the probability of assigning the spot i to cluster j based on the combined latent representation Y as follows: q i j = 1 + y i − μ j 2 / ρ − ρ + 1 2 ∑ j ′ 1 + y i − μ j ′ 2 / ρ − ρ + 1 2 (9) where μ j is the cluster center by K-means on learned representations. y i is the i th column of Y . We regard Q = q i j as the distribution of the assignments of all samples. ρ is the degree of freedom of student’s t-distribution.

To optimize the AE-GCN-learnable representation from the high-confidence assignment, we make data representation closer to cluster centers for improving the cluster cohesion. Hence, we calculate the target distribution P as follows: p i j = q i j 2 / s j ∑ j ′ q i j ′ 2 / s j ′ (10) where s j = ∑ i q i j is the soft cluster frequency. Each assignment in Q is squared and normalized to produce the target distribution P , which makes the data representation surround the cluster centers closer and helps AE-GCN learn a better representation for the clustering task. By minimizing the KL (Kullback–Leibler) divergence loss between Q and P distributions, the target distribution P can help the AE-GCN learn a better representation for the clustering task, i.e., making the data representation surround the cluster centers closer, thus leading to the following loss function: L c l = K L P Q = ∑ i ∑ j p i j log p i j q i j (11)

This design is regarded as a clustering-aware contrastive mechanism, where the P distribution supervises the updating of the distribution Q , and the target distribution P is calculated by the distribution Q in turn. Using this mechanism, AE-GCN can directly concentrate two different objectives: clustering objective and data reconstruction objective, in one loss function. Thus, the overall loss function of AE-GCN is L o b j = L r e s + β L c l (12) where β denotes the tunable parameter β > 0 and can be flexibly set, which balances data reconstruction and clustering optimization.

The top 3,000 highly variable genes (HVGs) for 13 10x Visium datasets, one ST dataset, and one Slide-seqV2 dataset are selected using scanpy.pp.highly_variable_genes() from the Scanpy Python package. The log-transformation of the expression profiles is performed using scanpy.pp.log1p() on the original gene expression data.

AE-GCN uses the combined latent representation Y to detect spatial domains by Leiden (Traag et al., 2019) algorithms implemented as scanpy.tl.leiden(). The parameter “resolution” can be adjusted to match the number of the manual annotations.

For the enhanced expression matrix X , AE-GCN aggregates the shared information between each spot and its surrounding neighbors by incorporating prior spatial information into gene expression, which is used to adjust expression values in each spot and enrich spatial local signals. For the reconstructed expression data X ′ , AE-GCN uses AE and GCN to reconstruct the enhanced expression matrix X . By minimizing the reconstruction error, the reconstructed data X ′ can reflect both the spatial local signals and expression measurement global signals. Thus, AE-GCN uses the reconstructed expression data X ′ as the denoised profiles.

We use adjusted Rand index (ARI) (Hubert and Arabie, 1985) and cluster purity (i.e., Eq. 13) (Zhao et al., 2021) to quantify the accuracy of the identified spatial domain and the reference annotations from original publications. c l u s t e r   p u r i t y = 1 N ∑ c ∈ C max g ∈ G c ∩ g (13) where C is denoted as the set of the spatial cluster set and G is regarded as the set of annotated groups. Due to cancer slices with rough annotations (e.g., IDC and PDAC cancer data), cluster purity is specifically used to evaluate the clustering performance on SRT cancer datasets (Zhao et al., 2021).

We use bulk expression data with patient survival information to evaluate the prognostic significance of genes via the Kaplan–Meier plotter (Zwyea et al., 2021) in the IDC and PDAC cancer studies.

We evaluated the ability of AE-GCN to detect spatial domains using 12 human dorsolateral prefrontal cortex (DLPFC) slices generated using 10x Visium. The DLPFC dataset obtained from spatialLIBD (Pardo et al., 2022) is manually annotated as the layered regions by gene markers and cytoarchitecture. The annotations can be considered as the ground truth for benchmarking. Based on this dataset, we compared AE-GCN with the existing state-of-the-art methods, including six spatial clustering methods [i.e., BayesSpace (Zhao et al., 2021), Giotto (Dries et al., 2021), SEDR (Fu et al., 2021), SpaGCN (Hu et al., 2021), stLearn (Pham et al., 2020), and STAGATE (Dong and Zhang, 2022)] and three non-spatial algorithms [i.e., variational autoencoder (VAE) (Kingma and Welling, 2019), Leiden implemented in Scanpy (Wolf et al., 2018), and Louvain implemented in Seurat (Butler et al., 2018)]. The adjusted Rand index (ARI) is used to quantify the similarity between the manual labels and identified clusters, which ranges from 0 for poor consistency to 1 for identical clusters.

Generally, most of the spatial clustering methods performed better than non-spatial algorithms (Wilcoxon signed-rank test P < 10 − 5 , Figure 2A), which showed that the integration of spatial information is necessary to improve the spatial clustering performance. Strikingly, AE-GCN had the highest mean ARI (mean ARI = 0.561) and substantially performed better than the competing methods over the slices (Wilcoxon signed-rank test P < 10 − 8 , Figure 2A). Taking slice 151673 as an example (Figure 2B), we found AE-GCN (ARI = 0.623), STAGATE (ARI = 0.588), BayesSpace (ARI = 0.556), and SEDR (ARI = 0.515) delineated the layered regions (Figure 2C). Notably, the partitions from AE-GCN (termed as a deep learning model-combined method) exhibited clearer and less noisy outcomes than those from single model-based methods (e.g., GCN-based SpaGCN and the VAE model).

To illustrate the effectiveness of AE-GCN on high-resolution SRT platforms, we applied AE-GCN to a mouse hippocampus Slide-seqV2 dataset (n = 41,786 spots). Slide-seqV2 can measure gene expression at near-cellular resolution (Stickels et al., 2021) but has lower number of transcripts per location/spot and higher dropouts than the 10x Visium platform. Thus, it poses more challenges for accurately distinguishing tissue structures from the data of high sparsity. To better validate the performance of AE-GCN, we also compared it with other domain detection methods and used the corresponding anatomical diagram from the Allen Mouse Brain Atlas (Sunkin et al., 2012) as the illustrative reference (Figure 3A).

Comparing with the reference, we found that AE-GCN and STAGATE can identify the spatially coherent domains compared to other involved methods. However, AE-GCN performed better to detect the fine-grained structures, such as the cornu ammonis 2 (CA2, AE-GCN domain 16), ventricle (AE-GCN domain 12), and habenula (AE-GCN domain 11) sections (Figure 3A). These sections are delineated with sharper boundaries and higher concordance with the anatomical annotation. We further isolated the focused regions and provided validations from other perspectives (Figures 3B–D). For the CA2 section, which is only detected by AE-GCN, the domain location showed good alignment with the marker gene expression (i.e., Pcp4 (San Antonio et al., 2014)) and independent in situ hybridization (ISH) image (Figure 3B). For ventricle and habenula sections, AE-GCN domains are closer to the shapes of their respective marker expression (Enpp2 for ventricle (Koike et al., 2006) and Gabbr2 for habenula (De Beaurepaire, 2018)) or stained regions and match the anatomical shape well (Figures 3C, D). Thus, for higher-resolution SRT data, AE-GCN is capable of effectively unveiling the fine-grained anatomical functional regions.

To illustrate the effectiveness of AE-GCN on cancer tissue, we applied AE-GCN to the human pancreatic ductal adenocarcinoma (PDAC) ST dataset (n = 428 spots). The histopathological image and annotations were taken as references (Figures 4A, B). We assessed these spatial domain identification methods using cluster purity (see Methods) as the quantitative measure on cancer datasets with rough annotation information. AE-GCN achieved the highest cluster purity (purity = 0.756) and detected more spatially enriched functional regions in tumor tissue than other compared methods (Figure 4C).

Next, we examined whether AE-GCN could provide more insights into the underlying tumor heterogeneity, as data sparsity could hinder other downstream analytical tasks, for example, the identification of differentially expressed genes (DEGs). In this manner, we used the AE-GCN-reconstructed data to denoise the low-quality measurements and evaluated the effectiveness in recovering gene spatial expression patterns. Based on the denoised data, we selected the top 50 DEGs of each domain from the reconstructed data X ′ and compared the log fold change (LFC) of these DEGs before and after denoising (Figure 4D). Overall, the comparison highlights the significant improvement of biological specificity brought by AE-GCN denoising across the identified domains (Wilcoxon signed-rank test P < 10 − 14 , Figure 4D). In particular, we found that some DEG expression (e.g., S100P and TNS4) appeared more spatially smoothed on spots in situ (Figure 4E). These two DEGs were validated to be the potential prognostic risk factors for PDAC (Figure 4F). For example, S100P is ever reported to be involved in the aggressive properties of cancer cells and associated with poor prognosis (Wang et al., 2012) (Figure 4F); TNS4 is associated with cancer cell motility and migration, whose high expression can indicate poor prognosis (Sakashita et al., 2008). These results indicate that AE-GCN has the potential to provide the in-depth biological insights into the underlying tumor heterogeneity from the perspectives of spatial domain detection and gene expression pattern recovery.

To illustrate the generalization ability of AE-GCN on cancer tissues, we next tested AE-GCN using the invasive ductal carcinoma (IDC) Visium dataset (n = 4,727 spots). The histopathological annotations from the original paper (Zhao et al., 2021) were taken as the reference (Figures 5A, B). We found that the identified domains of AE-GCN were highly consistent with the manual annotations (purity = 0.865, Figure 5C). Compared with the domains captured by other methods, the clustering partitions from AE-GCN showed clear spatial separations with few scatter points and high regional continuity.

Then, for functional gene identification, we identified the top 50 DEGs of each cluster from the denoised data X ′ . Similarly, based on the comparison before and after denoising, we found that AE-GCN significantly improves the LFCs of gene expression, revealing more biological specificity across domains, which may suggest the detection of new disease-associated genes (Figure 5D). For example, SLC7A5 and RDH16 are two newly found DEGs after denoising, whose spatial expression patterns are greatly enhanced after denoising (Figure 5E). Moreover, the two novel DEGs were shown to be the potential prognostic risk genes for breast cancer via survival analysis of independent clinical data (Figure 5F). Their biological functions in tumors indicate the prognostic relevance from previous studies. For example, SLC7A5 is reported to involve in tumor cell metabolism and promotes cell proliferation (El Ansari et al., 2018). RDH16 affects retinol metabolism to participate indirectly in breast cancer occurrence and progression (Gao et al., 2020). The application, along with the PDAC case, demonstrates that AE-GCN can unveil cancer heterogeneity from SRT data, enabling the discovery of novel spatial patterns of both samples and genes.