GraphPCA: a fast and interpretable dimension reduction algorithm for spatial transcriptomics data

Abstract

The rapid advancement of spatial transcriptomics technologies has revolutionized our understanding of cell heterogeneity and intricate spatial structures within tissues and organs. However, the high dimensionality and noise in spatial transcriptomic data present significant challenges for downstream data analyses. Here, we develop GraphPCA, an interpretable and quasi-linear dimension reduction algorithm that leverages the strengths of graphical regularization and principal component analysis. Comprehensive evaluations on simulated and multi-resolution spatial transcriptomic datasets generated from various platforms demonstrate the capacity of GraphPCA to enhance downstream analysis tasks including spatial domain detection, denoising, and trajectory inference compared to other state-of-the-art methods.

Background

Spatial transcriptomics (ST) technologies have fundamentally reshaped the current research landscapes of cellular and molecular biology, and significantly deepened our understanding on cellular heterogeneity, gene expression-cellular microenvironment interaction, and spatial specificity of gene expression in complex tissues [1,2,3,4]. Unlike single cell RNA sequencing (scRNA-seq) which loses spatial information during cell dissociation [5,6,7,8], ST quantifies expression of single RNA molecules while preserving spatial location information through in situ or spatial barcoding approaches [9]. Currently, there are two main types of ST techniques depending on the underlying experimental protocols, i.e., imaging-based methods and next-generation sequencing-based (NGS-based) technologies. The first type of methods includes in situ hybridization (e.g., seqFISH [10, 11], MERFISH [12,13,14], osmFISH [15]) and in situ sequencing (STARmap [16] and FISSEQ [17]). These methods often provide finer cellular positional details (in single-cell or even subcellular resolution), but are limited to pre-selected encoding probes, and thus fail to cover genes on the whole transcriptome scale [18]. The second type of techniques encodes location information into transcripts and quantifies both gene expression and spatial information through massively parallel sequencing, with representative technologies of ST [19], 10X Genomics Visium [20], Slide-seq [21], Slide-seqV2 [22], and Stereo-seq [23]. All these diverse technologies yield a data matrix that records the whole transcriptome expression levels at every spot, cell, or spatial pixel [3].

However, analyzing gene expression profiles obtained from various ST technologies faces enormous challenges due to the sparsity, ultra-high dimensionality, and low signal-to-noise ratio (SNR [24]) of the data matrix. Consequently, dimension reduction becomes a necessary preprocessing step to improve SNR and mitigate the curse-of-dimensionality [24, 25]. It is also a critical step in the analysis of scRNA-seq data, often followed by key downstream tasks such as cell type identification [26,27,28], trajectory inference [29,30,31], and batch effect removal [32,33,34,35]. With the rise of ST technologies, most researchers directly apply dimension reduction methods developed for bulk or scRNA-seq data to ST data (Seurat [36], Scanpy [37], STUtility [38]). However, those methods may fall short in terms of the capability of fully exploiting the location information in ST data, potentially leading to efficiency loss or even biased and erroneous results.

To overcome this limitation, dimension reduction methods crafted specifically for ST data have been proposed [39,40,41]. These methods typically assume that proximal spots exhibit similar gene expression patterns, such that the low-dimensional embedding can approximately preserve this spatial structure as much as possible. For instance, building upon the framework of probabilistic PCA, Shang and Zhou proposed SpatialPCA [39], which uses a kernel matrix to model the spatial correlation structure across locations. Also under the probabilistic PCA framework, DR-SC [40] integrates dimension reduction and clustering using a two-layer hierarchical Bayesian model, where the low-dimensional latent variables in it are assumed to follow a Gaussian mixture distribution. Both methods utilized computationally intensive iterative procedures, i.e., maximum likelihood-based optimization in SpatialPCA and MCMC in DR-SC, to approximate the locally optimal solution. There are also some recent approaches using deep learning algorithms to learn a nonlinear, low-dimensional embedding of the ST data. For example, SpaGCN [42] uses graph convolutional network to integrate multi-modal data including gene expression, spatial location of spots/cells, and histology data. STAGATE [43] obtains spatial representation of spots by a graph attention auto-encoder guided by a spatial neighbor network and cell type-aware network. However, deep learning-based methods suffer from lack of interpretability, high computational complexity in model training, and difficulty in tuning a large number of hyperparameters.

In this study, we develop GraphPCA, a novel graph-constrained, interpretable, and quasi-linear dimension-reduction algorithm tailored for ST data. GraphPCA learns the low-dimensional representation of ST data based on PCA with minimum reconstruction error, by incorporating spatial location information as constraints in the reconstruction step. By increasing the importance of the spatial network constraints, adjacent spots in the original dataset are more inclined to be positioned in nearby points in the low-dimensional embedding space. More importantly, the computationally efficient close-form solution of GraphPCA allows rapid embedding of massive single-cell or even subcellular resolution spatial transcriptomic data generated from techniques including Slide-seq and Stereo-seq. Finally, we demonstrate the superiority of GraphPCA over competing methods, i.e., PCA, NMF, SpaGCN, BayesSpace [44], DR-SC, SpatialPCA, and STAGATE, through comprehensive synthetic experiments and real spatial transcriptomic data with a variety of resolutions, species, and tissue states.

Results

Overview of GraphPCA

Here, we developed GraphPCA, an interpretable and quasi-linear statistical algorithm for dimension reduction of spatial transcriptomics data (Fig. 1a, see “Methods” for details). Building upon the flexible PCA framework, GraphPCA enables the low-dimensional embeddings to effectively preserve location information by leveraging spatial neighborhood structure between spots/cells as graph constraints. The input data for GraphPCA includes a gene expression matrix along with spatial coordinates of spots, which are utilized to construct a spatial neighborhood graph ($k$ NN graph by default [45]). In contrast to the classical PCA, GraphPCA infers an embedding matrix integrating both spatial location and gene expression information by solving an optimization problem with constraints determined by the constructed spatial neighborhood graph. In particular, we incorporated the constraints by introducing a penalty term (Fig. 1a, middle panel), with penalty strength controlled by a tunable hyperparameter λ. By design, the resulted constrained penalized optimization problem has a closed-form solution. As a consequence, GraphPCA can process ST data efficiently at vastly different scales. The low-dimensional spatial embeddings inferred by GraphPCA can be readily utilized for various downstream analysis tasks including spatial domain detection, visualization, denoising, and trajectory inference (Fig. 1a, right panel).

Fig. 1

Workflow of GraphPCA and synthetic experiments validation. a GraphPCA is a novel graph-constrained, interpretable, and quasi-linear dimension-reduction method tailored for ST data. It begins with input data including gene expression matrix X along with spatial coordinates S of spots. It first constructs a spatial neighborhood graph G using the spatial coordinates S and sets the graph constraint parameter $\lambda$ to characterize the spatial relationships and dependences in the low-dimensional embedding Z. Subsequently, GraphPCA infers the embedding matrix Z by integrating both spatial location S and gene expression information X by solving a non-convex optimization problem with graph constraints. The output of GraphPCA can be readily utilized for various downstream analysis tasks including spatial domain detection, trajectory inference, and denoising. b In simulation, we obtained the anatomical structure of mouse brain sagittal from the Allen Brain Atlas as ground truth layer labels and simulated ST data using scDesign3. c Clustering accuracies of GraphPCA on simulated datasets across varying values of the graph constraint parameter $\lambda$ (x-axis). d Clustering accuracies of different methods on simulated data. Error bars indicate 95% confidence intervals across 20 replicates. e The robustness of GraphPCA and other methods under varying simulation scenarios including different sequencing depths, noise levels, spot sparsity, and expression dropout rates. For each scenario and method, dots represent the mean ARI calculated across 20 replicates

Synthetic experiments

To comprehensively evaluate the performance of GraphPCA and other competing algorithms, we generated a series of simulated data using scDesign3 [46]. With simulated datasets, one could use curated spatial domain labels as ground truth to evaluate the accuracy of the clustering results. Specifically, we downloaded the anatomical structure of mouse brain sagittal from the Allen Brain Atlas [47] as ground truth layer labels (Fig. 1b) and used scDesign3 to simulate gene expression data (see “Methods” for details). We compared GraphPCA with four popular algorithms for ST dimension reduction (PCA, NMF, DR-SC, and SpatialPCA) and three algorithms for spatial domain detection (SpaGCN, STAGATE, and BayesSpace). Clustering performance was evaluated using the adjusted Rand index (ARI [48]), normalized mutual information (NMI [49]), and homogeneity score (HS [50]).

We first examined the impact of the hyperparameter $\lambda$ on spatial domain detection. By increasing $\lambda$ from 0 to 1 with a step size of 0.01 across four synthetic datasets, we first observed an overall improvement in ARI as $\lambda$ increases (Fig. 1c). This suggests that the integration of proper location information enhances spatial domain detection, with a larger graph constraint resulting in smoother and more contiguous clustering output (Additional file 1: Figs. S1–2). However, excessively large $\lambda$ leads to the spatial constraints dominating the dimension reduction objective and results in deteriorated performance. Based on the empirical observations from the synthetic experiments, we recommend choosing $\lambda$ between 0.2 and 0.8 (default 0.3 unless stated otherwise) for tissue samples with evident spatial layered structure. Compared to other methods, GraphPCA demonstrates superior performance on the synthetic data (median ARI: 0.784), outperforming algorithms considering spatial information such as PCA (median ARI: 0.556) and NMF (median ARI: 0.185) across all metrics (Fig. 1d, Additional file 1: Fig. S3a). Methods incorporating spatial information (DR-SC and SpatialPCA) and deep learning-based algorithms (SpaGCN and STAGATE) also achieve comparable clustering results, but exhibit relatively high variabilities (standard deviations across 4 replicates are 0.040, 0.053, 0.088, 0.071, and 0.044 for GraphPCA, DR-SC, SpatialPCA, SpaGCN, and STAGATE, respectively).

Next, we evaluated the robustness of GraphPCA under varying simulation scenarios including different sequencing depths, noise levels, spot sparsity, and expression dropout rates (Fig. 1e, Additional file 1: Fig. S3b). For sequencing depths, we down-sampled the simulated count matrix from 100 to 10% of the original depth and found that GraphPCA consistently outperformed other competing methods even with only 10% of the original sequencing depth (Fig. 1e, first column). By introducing different levels of Gaussian white noises with increased standard deviations to the raw data, GraphPCA performs consistently well despite the increasing SNR, demonstrating the robustness of GraphPCA to noise. In contrast, ST-tailored algorithms such as SpaGCN, DR-SC, and BayesSpace eventually fail to cluster cells correctly under low SNRs (Fig. 1e, second column). For spot sparsity, we randomly removed a fraction of spots to increase sparsity of the data (from 10 to 90%) and found that both GraphPCA and BayesSpace sustain relatively high performance under all sparsity ratios (Fig. 1e, third column). To examine the impact of expression dropout, we randomly set the expressions of a fraction of genes (from 10 to 90% of all genes) to 0 and found that GraphPCA still maintains high ARIs even at a dropout rate of 60% (Fig. 1e, fourth column). In summary, GraphPCA exhibits superior accuracy and robustness in clustering performance across all four scenarios compared to competing methods.

Human dorsolateral prefrontal cortex data by Visium

We next assessed the clustering performance of GraphPCA against competing methods using real ST data with expert pathological annotations. These datasets encompass various sequencing technologies, resolutions, species, and tissue states. We first evaluated GraphPCA on the human dorsolateral prefrontal cortex data (DLPFC) obtained from the 10X Visium platform [51]. This dataset comprises twelve tissue sections from three adult donors, with a median depth of 291 million reads per sample. Each tissue section has a median of 3844 spots (Additional file 1: Table S1), and on average 1734 genes are recorded per spot. As an example, the slice 151673 assays 3639 spots across 33,538 genes, and each spot is annotated into one of six neuronal layers or white matter. We compared GraphPCA with other methods, using the true number of spatial domains and the same 3000 spatially variable genes (SVGs) identified by SPARK [52] as input.

As expected, dimension reduction algorithms that do not consider spatial information, i.e., PCA and NMF, fail to distinguish distinct layers of cortex, thus resulting in poor clustering accuracies (PCA: ARI 0.290, NMF: ARI 0.235) (Fig. 2a). Although incorporating additional histological imaging data as input, SpaGCN fails to accurately recover cortical layers 3–6. BayesSpace does not capture the boundary between white matter and layer 6 correctly. DR-SC identifies cortical layers that roughly match manual annotations in shape but display discrepancies in thickness, and a few spots adjacent to white matter are erroneously identified as layer 1 (Fig. 2a). In contrast, GraphPCA, STAGATE, and SpatialPCA effectively delineate cortical layers and accurately recover layer boundaries. GraphPCA achieves superior performance among compared methods, as indicated by higher ARI, NMI, and HS scores (0.536, 0.689, and 0.706, respectively), which nearly double those by space-unaware methods (PCA: ARI 0.290, NMI 0.452, HS 0.456; NMF: ARI 0.235, NMI 0.292, HS 0.287). Sensitivity analysis indicates that GraphPCA obtains stable spatial domain detection accuracy over varying numbers of SVGs for DLPFC data (Additional file 1: Fig. S4). While STAGATE (ARI 0.576) and SpatialPCA (ARI 0.561) obtain slightly higher ARI, they are significantly more time-consuming due to employing the iterative optimization procedures (Fig. 2c). More importantly, GraphPCA has the distinct advantage of providing an analytic solution, making it more interpretable and robust compared to deep learning-based algorithms like SpaGCN and STAGATE. Similar results are also observed in other sections of the DLPFC data (Fig. 2b, Additional file 1: Fig. S5).

Fig. 2

Performance of GraphPCA on human dorsolateral prefrontal cortex (DLPFC) samples. a Clustering results of various methods on sample 151673 of the DLPFC dataset. The ground truth of tissue regions is the manual annotation of six cortex layers and white matter (WM), as provided by the original study. Manual annotations and clustering results of other DLPFC slices are shown in Additional file 1: Fig. S5. b Boxplots illustrating clustering accuracies across all 12 sections by eight methods. c Runtimes of different methods for spatial domain detection across all 12 sections. d UMAP visualizations of the sample 151673 generated by GraphPCA, PCA, NMF, and SpatialPCA embeddings, respectively. e Line plots showing the dynamic change in clustering accuracy at varying numbers of input low-dimensional components (x-axis) inferred by GraphPCA and other methods. f Scatter plot displaying the mean expression of layer-specific markers including MBP and CLDND1 (white matter), KRT17 (layer 6), PCP4 (layer 5), PVPLB and TMSB10 (layer 4), MFGE8 (layer 3), CXCL14 (layer 2), and MALAT1 and FABP3 (layer 1). Clusters in y-axis correspond to the spatial domains inferred by GraphPCA. g Trajectory inference on sample 151673 based on the inferred low-dimensional components of GraphPCA. Arrows point from tissue locations with low pseudo-time to those with high pseudo-time

We then performed UMAP visualization of low-dimensional embeddings generated by each method. It is shown that different spatial domains inferred by GraphPCA are clearly separated, with their relative orders and global arrangements aligned well with histological images (Fig. 2d, Additional file 1: Fig. S6). This observation suggests that the GraphPCA approach effectively captures the local structure and relative positioning of the cells in the embedding space, which helps provide valuable insights into the spatial organization of the tissue. In contrast, spots associated with different layers are intermixed in UMAP by methods without explicitly considering spatial constraints (Fig. 2d, Additional file 1: Fig. S6). We next examined the impact of the number of principle components (PCs) on the clustering accuracy. Since GraphPCA provides a globally optimal solution without relying on iterative procedures and random initial values, it exhibits more stable clustering results across different numbers of PCs in terms of both ARI and McFadden-adjusted pseudo-R² (Fig. 2e, Additional file 1: Fig. S7). In contrast, NMF, SpaGCN, STAGATE, DR-SC, and SpatialPCA exhibit larger variabilities in clustering accuracies, indicating their high sensitivity to the embedding dimension and parameter initialization. Additionally, the seven spatial domains identified by GraphPCA are enriched with layer-specific marker genes reported previously [51] (Fig. 2f, Additional file 1: Fig. S8). We also performed gene set enrichment analysis (GSEA) to identify significant pathways that are highly enriched in each layer based on these layer-specific genes. The enriched pathways include neuron projection, cytoplasmic translation, and synapse-related, which represent common and crucial transcriptional programs underlying cortical development and maturation (Additional file 1: Fig. S9). Further analysis on the top PCs by GraphPCA reveals associations with spatial expression patterns of white matter, layer 3, and layer 4 (Additional file 1: Fig. S10a). The top 5 weighted genes by absolute values in each PC exhibit consistent spatial patterns, validating that the extracted patterns represent co-expression modules in corresponding regions (Additional file 1: Fig. S10b–d). Trajectory analysis by GraphPCA identifies an inward-outward spatial trajectory from white matter to cortical layers (Fig. 2g, Additional file 1: Fig. S11a), which agrees with previous studies that new neurons are generated in the ventricular zone and migrate outwards along radial glial fibers to integrate into existing layers [53, 54]. In contrast, the trajectory patterns produced by PCA, SpaGCN, STAGATE, and DR-SC appear to be random, failing to capture coordinated spatial-gene expression relationships (Additional file 1: Fig. S11b, c).

Mouse medial prefrontal cortex data by STARmap

Next, we applied GraphPCA to the high-resolution image-based ST data, specifically the mouse medial prefrontal cortex (mPFC) data generated by STARmap [16]. This dataset consists of 1049 cells measured across 166 genes, providing single-cell resolution for each spot. The mPFC is a vital cognitive control region that potentially involved in anxiety and fear regulation in human and mouse, and previous studies have divided into four laminar domains, i.e., L1, L2/3, L5, and L6 (Fig. 3a). We downloaded cell-type annotation of each spot from the original paper [16] and domain-level annotation in [55] as ground truth to evaluate clustering performance (Fig. 3b, g).

Fig. 3

Evaluation of GraphPCA on mouse medial prefrontal cortex (mPFC) data generated by STARmap. a An anatomic reference atlas displays the structure of the prelimbic area in the mouse prefrontal cortex. b Manual annotation provided by the original study. c Spatial domains detected by GraphPCA. d Barplots illustrating spatial domain detection ARI by different methods on the mPFC data. e Visualization of spatial domains identified by PCA, NMF, DR-SC, and SpatialPCA. The clustering results of SpaGCN, BayesSpace, and STAGATE are shown in Additional file 1: Fig. S12a. f Trajectories of the mPFC tissue based on low-dimensional embeddings obtained by GraphPCA, PCA, SpaGCN, and SpatialPCA. Trajectories are inferred by Slingshot under default parameters. g Annotated cell type labels for the mPFC data are available from the original study. h Barplots illustrating cell type clustering ARI by different methods on the mPFC data. i Spatial composition of cell types identified by different methods. Excitatory neurons: eL2/3, L5-1, eL5-2, eL5-3, eL6-1, and eL6-2; inhibitory neurons: Reln, VIP, SST, NPY, and Lhx6; Oligo: oligodendrocytes; Smc: smooth muscle cells; Astro: astrocytes; Endo: endothelia cells; Mixed: group of cells expressing marker genes associated with multiple cell types; and Unidentified: cells with no clear marker gene expression. j Cell type identification results by GraphPCA and other methods. Top: barplots displaying the number of cell types identified by each method. Bottom: heatmap showing correlations between estimated cell type proportions from each method and the ground truth

We first performed dimension reduction on the expression data and examined their ability to recover the spatial domain structure. As anticipated, GraphPCA accurately identifies four spatial domains, closely matching the ground truth and achieving the highest clustering scores (ARI: 0.886; NMI: 0.85; HS: 0.848) (Fig. 3c, d, Additional file 1: Fig. S12). The high performance of our method is consistent across different numbers of input low-dimensional components, as evidenced by ARI and McFadden-adjusted pseudo-R² metrics (Additional file 1: Fig. S12d). In contrast, clustering results reveal that PCA, NMF, SpaGCN, BayesSpace, and SpatialPCA struggle to effectively distinguish cortical layers L2/3 and L5, while incorrectly identifying layer L1 as a discrete structure in space. Although DR-SC yields relatively compact clusters, it fails to identify L2/3 (Fig. 3e). STAGATE misidentified a few spots of layer L6 as layers L1 and L5 (Additional file 1: Fig. S12a). The precise layer detected by GraphPCA allows us to identify layer-specific marker genes. We then performed differential gene expression analysis and revealed known layer marker genes including Bgn (L1), Cux2 (L2/3), Tcerg1l (L5), and Pcp4 (L6) (Additional file 1: Fig. S13). Furthermore, GraphPCA takes the least runtime, even surpassing PCA, making it highly scalable for the incoming large-scale data (Additional file 1: Fig. S12c). Finally, trajectory inference using top PCs obtained from different dimension reduction methods indicates that GraphPCA correctly identifies a path from L1 to L6, aligned well with known cortical development, while other methods fail to do so (Fig. 3f, Additional file 1: Fig. S14).

We next performed cell type clustering on the same dataset, in which every cell is annotated as one of the 15 types including excitatory neurons: eL2/3, L5-1, eL5-2, eL5-3, eL6-1, and eL6-2; inhibitory neurons: Reln, VIP, SST, NPY, and Lhx6; Oligo: oligodendrocytes; Smc: smooth muscle cells; Astro: astrocytes; and Endo: endothelia cells (Fig. 3g). We compared GraphPCA with two popular cell type clustering methods for scRNA data (Seurat [36] and SC3 [56]) and two methods designed for single-cell resolution ST data (BASS [55] and FICT [57]). To obtain the clustering results in cellular level, we reduced the graph regularization strength of GraphPCA to increase the clustering resolution ($\lambda$ = 0.2, parameter settings for different datasets and analysis tasks are listed as Additional file 1: Table S2). GraphPCA achieves the best clustering performance (ARI: 0.436), surpassing other methods (Seurat: 0.29; SC3: 0.225; BASS: 0.414; and FICT: 0.245) when the true number of cell types is specified in the subsequent K-means clustering (Fig. 3h). Similar conclusions can be drawn when evaluated by NMI and HS (Additional file 1: Fig. S15b).

Then, we analyzed the spatial distribution of cell types identified by various methods (Fig. 3i). GraphPCA successfully detects the known enriched cell types in their corresponding layers, such as Smc and Endo cells in L1, eL2/3 cells in L2/3, eL5-3 and L5-1 cells in L5, and eL6-1 and eL6-2 cells in L6. It accurately distinguishes the excitatory neurons at layer 6 (eL6) into two subtypes (eL6a and eL6b) and identifies inhibitory neurons VIP and its subtype Lhx6 with clear marker gene expression. Overall, GraphPCA identifies the highest number of correct cell types, while other methods produce multiple unidentified cell clusters without known marker genes, or mixed cell clusters with ambiguous marker gene expressions (Fig. 3j, Additional file 1: Figs. S16–23). In summary, the embeddings derived by GraphPCA effectively capture cell-identity differences at both domain and cell-type levels (Additional file 1: Figs. S24, 25).

We next extended the application to other NGS-based ST technologies with subcellular resolution. Based on the Slide-seq data of mouse cerebellum [21] and Slide-seqV2 data of mouse hippocampus [22], GraphPCA demonstrates superior clarity and continuity in representing the Purkinje layer in the cerebellum data (Additional file 1: Fig. S26). For the hippocampus data comprising expression of 23,264 genes in 53,208 spots, GraphPCA clearly represents the hippocampus, dentate gyrus, and CA1 and CA3 layers, which are not clearly defined by other methods (Additional file 1: Fig. S27).

Murine liver data by 10X Visium

To validate the performance of GraphPCA on more complex tissues, we next applied it to the murine liver data generated by 10X Visium [58], which consists of 1293 cells measured across 31,053 genes. The murine liver is a heterogeneous tissue composed of hexagonal lobules that are radially polarized by blood flow and morphogens [59, 60]. As blood flows directionally toward the central vein, hepatocytes take up oxygen and nutrients, metabolize hormones, and create a gradient along the periportal-pericentral axis known as “liver zonation” [60,61,62] (Fig. 4a). According to the zonation annotation by Guilliams et al. [58], the murine liver lobule can be divided into four zones, i.e., central vein, mid zone, periportal zone, and portal vein (Fig. 4b).

Fig. 4

Evaluation of GraphPCA on murine liver data generated by 10X Visium. a An anatomic reference atlas displays the structure of the murine liver lobule area. b Manual annotation provided by the original study. c Spatial domains detected by GraphPCA. d Visualization of spatial domains identified by SpaGCN, STAGATE, and SpatialPCA. The clustering results of PCA, NMF, BayesSpace, and DR-SC are shown in Additional file 1: Fig. S28a. e Barplots illustrating spatial domain detection ARI by different methods on the murine liver data. f Expression of zonation marker genes identified by GraphPCA in the ROI. g Line plots illustrating expression level of zonation marker genes identified by GraphPCA along the portal–central lobule axis (x-axis) in the ROI. Error bars indicate 95% confidence intervals across each layer spots. h Spatial trajectory and pseudo-time of the ROI based on low-dimensional embeddings obtained by GraphPCA. Trajectories are inferred by Slingshot under default parameters. Arrows point from tissue locations with low pseudo-time to those with high pseudo-time. Color represents different tissue regions

We first performed dimension reduction and spatial domain detection using different methods on the murine liver data. GraphPCA accurately identifies each zone along the lobular axis in the correct sequential order: portal vein, periportal zone, mid zone, and central vein (Fig. 4c). However, other space-aware methods struggle to identify the portal vein within the periportal zone (Fig. 4d, Additional file 1: Fig. S28a). For instance, STAGATE and SpatialPCA misclassify the portal zone as part of the periportal zone, while BayesSpace produces overly smoothed results, leading to the emergence of two dominant clusters and failing to recognize the characteristic zones of the hepatic lobule. Overall, GraphPCA achieves the highest clustering performance (ARI: 0.622), which is at least 10% higher than those by other methods (PCA: 0.504; NMF: 0.027; SpaGCN: 0.413; STAGATE: 0.382; BayesSpace: 0.204; DR-SC: 0.22; and SpatialPCA: 0.207) (Fig. 4e).

Then, we selected the distinct zonal layers of the hepatic lobule in the tissue regions that were accurately detected by GraphPCA and performed differential expression analysis to identify zonation marker genes (Fig. 4f, Additional file 1: Fig. S29). The identified layers by GraphPCA are enriched with known marker genes (portal vein: Gpx3, Tmsb10, Gsn; central vein: Cyp2e1, Mup17, Mgst1) [59]. In particular, Gpx3 is an antioxidant enzyme that plays a crucial role in protecting cells from oxidative stress by reducing hydrogen peroxide and lipid peroxides [63,64,65]. Cyp2e1 is a member of the cytochrome P450 family involved in xenobiotic metabolism [66,67,68]. We noticed that Cyp2e1, Mup17, and Mgst1 are not detectable near annotated portal veins, whereas Cyp2f2, Alb, and Spp1 show high expression levels near the portal vein, with no signals in the central vein. Additionally, we visualized the spatial expression of zonation marker genes identified by GraphPCA for each zone (Fig. 4g), revealing a pronounced expression gradient along the lobular axis, consistent with a previous study [59].

We next performed GSEA to identify significant pathways that are significantly enriched in each layer based on these layer-specific genes (Additional file 1: Fig. S30). Pathway analysis of the portal vein marker genes reveals the strongest enrichment in genes associated with the process of “collagen-containing extracellular matrix.” This pathway not only provides essential physical scaffolding for cellular constituents but also initiates crucial biochemical and biomechanical cues necessary for tissue morphogenesis, differentiation, and homeostasis [69, 70]. Pathways related to antigen processing and presentation, specifically the “MHC class II protein complex” and “multivesicular body,” are also highly enriched within the portal vein, highlighting their role in influencing immune recognition and response [71]. In contrast, pathways “peroxisomes” and “microbodies” are associated with the central vein, underscoring their importance in maintaining the metabolic functions of hepatic tissue [72].

Finally, we performed trajectory inference by using top PCs obtained from different dimension reduction methods to identify one trajectory per method from portal vein to central vein (Fig. 4h, Additional file 1: Fig. S31). Specifically, GraphPCA captures the well-known spatial pattern of hepatocyte differentiation along the lobular axis [59, 73, 74]. As hepatocytes differentiate from the portal vein, they are gradually exposed to higher levels of oxygen, nutrients, and metabolites, creating a differentiation gradient that reflects the physiological maturation of hepatocytes in response to blood flow [60, 61, 74,75,76]. In contrast, the trajectories inferred from other space-aware methods appear to be random, with pseudo-time values intermingled across layers (Additional file 1: Fig. S31).

Taken together, the low-dimensional representation obtained from GraphPCA is accurate and robust in analyzing heterogeneous murine liver sample, enabling the exploration of transcriptional and functional heterogeneity across zones along the lobular axis between the portal and central veins.

Denoising gene expression profiles for deciphering gene spatial patterns

Leveraging the interpretability inherent in the PCA-based framework, GraphPCA is also capable of reconstructing the gene expression matrix from low-dimensional embeddings and projection vectors, thereby achieving denoising of the original expression data (see “Methods” for details). In this part, we aimed to unveil intrinsic spatial gene expression patterns in denoised DLPFC data using GraphPCA.

We conducted a comparative analysis involving GraphPCA and Sprod, state-of-the-art denoising methods for ST data. Specifically, we examined layer-marker genes of six cortical layers in sample 151673, as selected by Maynard et al. [51], and analyzed their spatial expression patterns before and after denoising by two methods (Fig. 5a). For both methods, the spatial expression of each layer-marker gene is revealed, with high expression regions matching the corresponding layers compared with raw data. Notably, in GraphPCA-denoised data, LAMP5 and NTNG2, marker genes for temporal and visual cortices, respectively, exhibit differential expression between layers 5 and 6, which is obscured in the raw data. While Sprod also recovers these spatial patterns, gene expression signals are subtle, lacking sufficient differential expression between layers. We further assessed the quality of the denoised data by two methods, as well as the raw data, measured by quantitative metrics including log fold change and log adjusted p values of known layer marker genes. GraphPCA achieves better performance compared to Sprod as indicated by higher log fold changes and negative log adjusted p values (Additional file 1: Fig. S32). In summary, GraphPCA improves the detection of spatial patterns and inter-layer expression differences of layer-marker genes (Fig. 5b), demonstrating its utility in denoising and elucidating intrinsic spatial gene expression patterns.

Fig. 5

GraphPCA deciphers gene spatial patterns of layer-specific marker genes in the DLPFC tissue. a Visualizations of raw and denoised spatial expressions of six layer-specific marker genes on sample 151673 of the DLPFC dataset. Benchmark methods include GraphPCA and state-of-the-art denoising tool Sprod. b Gene expression heatmaps of layer-specific marker genes on raw and denoised data generated by different methods

Integrating multiple ST datasets improves the performance of downstream analyses

As illustrated in previous sections, analysis based on single ST datasets can be inaccurate due to the high level of noise. In addition to denoising, another alternative strategy to enhance the signal is integrating multiple slices from the same sample or tissue type [55, 77]. To enable this function, we further designed GraphPCA_multi, a multi-sample extension of GraphPCA through aggregating multiple spatial neighborhood graphs and expression profiles into one matrix. The low-dimensional embedding per sample can be derived by optimizing the overall objective function across all samples (see “Methods” for details).

To validate the capability of our multi-sample integration model, we applied GraphPCA_multi to aggregate slice 151674 with three adjacent sections (slices 151673, 151675, 151676) and re-detected spatial domains using the joint embeddings. Clustering performance for slice 151674 significantly improves (ARI increased from 0.38 to 0.54), and the deteriorated pattern of layers 3 and 4 in single-slice estimation is recovered (Fig. 6a). Taking all samples, GraphPCA_multi outperforms single-sample analysis in terms of ARI (Fig. 6b) and surpasses state-of-the-art ST data integrating tools, BASSMult (a multi-sample version of BASS [55]) and STAligner [77] for most samples (Fig. 6c, Additional file 1: Fig. S33). In addition, GraphPCA_multi also achieves the lowest running time while exhibiting lower variance across DLPFC slices compared to BASSMult.

Fig. 6

Performance of GraphPCA in integration of multiple ST datasets. a Comparison of identified spatial domains by GraphPCA and GraphPCA_multi (the multi-sample version of GraphPCA) for sample 151674 of the DLPFC dataset. b Boxplot showing clustering accuracy of all 12 sections by GraphPCA, GraphPCA_multi, and two other multi-sample integration methods, i.e., STAligner and BASSMult. All three multi-sample methods take four neighboring sections for each section as input. Two stars (**) denote that the p value is less than 0.01 by paired one-sided t-test. c Inferred spatial domains on four slices (151507, 151508, 151509, and 151510) of DLPFC by GraphPCA_multi, STAligner, and BASSMult, respectively. d Inferred spatial domains on three sections (BZ5, BZ9, and BZ14) of the mPFC dataset by GraphPCA_multi, STAligner, and BASSMult, respectively. Pathological annotation by the original STARmap study is served as the ground truth. e UMAP visualization of low-dimensional embedding by GraphPCA_multi colored by three sections. f UMAP visualizations of low-dimensional embedding by GraphPCA_multi (top) and STAligner (bottom) colored by spatial domains

Next, we examined the performance of GraphPCA_multi in integrating samples from the same tissue type, not necessarily adjacent slides or from the same sample. To achieve this, besides the BZ5 sample used in pervious analyses, we obtained two additional tissue sections BZ9 and BZ14 generated by STARmap from different mice. Cells in all sections are annotated into four cortical layers (L1, L2/3, L5, and L6) according to the original publication [16]. As observed, GraphPCA_multi achieves superior spatial domain identification versus BASSMult and STAligner across all three sections. In detail, the four layers identified by GraphPCA_multi closely match the manual annotation in terms of tissue boundaries and thickness (Fig. 6d). In UMAP visualization using the top PCs from GraphPCA_multi, three tissue sections are evenly mixed (Fig. 6e), while different layers are clearly separated (Fig. 6f), indicating its potential in batch effect removal. Notably, this segregation accords with the functional similarity of adjacent layers and developmental trajectories. Similar to the DLPFC, the captured trajectory pattern reflected the developmental process of cortical layers, transitioning from existing layers to newer ones (from L1 to L6). In contrast, STAligner arranges cells in concentric clusters rather than matching organ developmental trajectory (Fig. 6f). These results validate the efficiency of GraphPCA_multi in integrating multiple samples to enable more accurate biological discovery.

Discussion

In this study, we presented GraphPCA, a novel dimension-reduction method tailored for ST data by combining the strengths of graphical constraint and PCA-based framework. By leveraging the graphical constraint, GraphPCA ensures that projections of adjacent spots are closer in the low-dimensional space. Each embedding dimension is highly associated with specific spatial gene expression patterns, enabling the gene-component projection matrix to reflect the differential spatial distribution of co-expressed gene modules. Through extensive simulations and validations on real datasets, we demonstrated that the low-dimensional embeddings by GraphPCA can enhance performance of downstream analyses including spatial domain detection, trajectory inference, and denoising. The model flexibility of GraphPCA allows us to easily extend it for multi-sample integration, which further improves clustering accuracy by incorporating gene expression information from other slices. Importantly, the closed-form solution exploited in GraphPCA ensures great computational efficiency compared to deep learning-based approaches, madding it easily scalable to the escalating throughput of emerging spatial transcriptomic assays (e.g., Stereo-seq, Slide-seqV2) that involve tens of thousands of spots.

The hyperparameter $\lambda$ in GraphPCA balances the relative importance of reconstruction error in dimensionality reduction and the spatial constraints, and thus a proper value of $\lambda$ is critical. When $\lambda$ = 0, GraphPCA reduces to the typical PCA; as $\lambda$ increases, GraphPCA embeddings incorporate more spatial information, but excessively large $\lambda$ may lead to less informative embeddings dominated by the spatial constraints. Therefore, we thus suggest that users try different values of $\lambda$ to evaluate the consistency of the clustering results, and carefully check if the results are biologically plausible with background knowledge in downstream data analyses. In practical applications, it is also helpful to generate multiple embeddings across a range of $\lambda$ and integrate all outcomes to improve robustness for downstream analyses. Other than spatial locations, the spot similarity network based on gene expression or H&E image has not been incorporated. In theory, we could impose different $\lambda$ to adjust for other spatial-related constraints, in such case, our model still has a closed-form solution (see “Methods” for details). The preliminary result indicates that integrating both spatial location and gene expression data provides a modest clustering performance boost over the spatial position-only model on the mPFC data (Additional file 1: Fig. S34).

The construction of spot neighborhood graph is another pivotal factor determining the quality of the low-dimensional embedding. We primarily considered $k$-nearest neighbor ($k$ NN) graphs based on the Euclidean distance, while the GraphPCA package also provides various other metrics to construct spatial neighborhood graph of spots/cells based on spatial coordinate information (e.g., circular-neighborhood graph (CG), Delaunay triangulation graph (DTG) [78], and shared nearest neighbor graph (SNN) [79], see “Methods” for details). For 10X Genomics Visium and Stereo-seq that exhibit grid-like spot arrangement patterns, these spatial graph construction methods can be made equivalent when setting tuning parameters to certain values. But for single-cell resolution ST techniques such as STARmap and MERFISH, different methods may produce different spatial neighborhood graphs and thus different low-dimensional embeddings. We also tested different spatial graph construction methods as well as their corresponding hyperparameters on the performance of our method based on 10X Genomics Visium and single-cell resolution STARmap data (Additional file 1: Fig. S35). We could observe a clear improvement in clustering accuracies as the connectivity of spatial graphs increased (more neighbors in $k$ NN and SNN and larger radius in circular-neighborhood graph), which suggests the important contribution of graphical constraints in enabling accurate spatial domain detection. However, excessively large spatial graphs can also be detrimental, as they may produce over-smoothed domain detection results (Additional file 1: Fig. S36).

Conclusions

In conclusion, we proposed an interpretable and quasi-linear dimension reduction algorithm, GraphPCA, that efficiently captures the biological signal from ST data. This scalable and extensible framework facilitates the learning of low-dimensional embeddings while integrating gene expression with spatial location, thereby enabling various downstream analyses. Through combining graph representation learning and the PCA framework, GraphPCA retains the advantage of linear embedding while capturing nonlinear spatial information underlying the tissue architecture. In the future, we plan to explore the adaptability of GraphPCA to other spatial techniques, such as the emergent spatially resolved proteomics, epigenomic, and metabolomic data, with the potential to unveil novel biological insights into the spatial orchestration of these crucial omics layers.

Methods

Data preprocessing

GraphPCA takes the gene expression matrix and spatial coordinate information from spatial transcriptomics data as input. Let $X$ be an $n\times m$ gene expression matrix, where ${x}_{ij}$ denotes the expression of the $j$ th gene on the $i$ th spot, $i=1,\cdots , n$, $j=1,\cdots , m$. The spatial coordinates of individual spots are denoted as $S={\left({s}_{1},\dots , {s}_{n}\right)}^{\intercal }, {s}_{i}\in {R}^{2}/{R}^{3}$. Before performing GraphPCA, we assumed that the gene expression counts have already been preprocessed with analytic Pearson residuals proposed by Lause et al. [80] and further scaled for each gene to have zero mean and unit standard deviation. Then, we selected the top 3000 spatial variable genes using SPARK package.

Construction of spatial neighborhood graph

Under the assumption that spots in close spatial proximity often exhibit similar gene expression patterns, we constructed a sparse spatial neighborhood graph to capture the complex spatial interactions based on the spatial coordinate information $S$. By default, an undirected graph is built by employing the $k$-nearest neighbor ($k$NN) algorithm based on Euclidean distance, where $k$ is the number of neighbors. Then, a sparse adjacent matrix $A$ based on $k$NN graph can be calculated as:

$${A}_{ij}=\left\{\begin{array}{c}1, j\in N\left(i\right),\\ 0,otherwise,\end{array}\right.$$

where $N(i)$ represents the set of neighbors of spot $i$. To characterize distinct spatial distributions of spots across datasets, GraphPCA also accommodates other spatial neighborhood graph construction approaches such as SNN graph, circular-neighborhood graph, and Delaunay triangulation [81].

Simulation design

We conducted a series of simulations to comprehensively evaluate the performance of GraphPCA and other dimension reduction and spatial domain detection algorithms. To achieve this, we downloaded mouse brain sagittal datasets from 10X Genomics and obtained the anatomical structure of the mouse brain sagittal from the Allen Brain Atlas as ground truth layer labels.

We utilized scDesign3, a multi-omics data simulator based on the generalized additive model for location, scale, and shape (GAMLSS), to generate the simulated gene expression matrix. scDesign3 models the joint distribution of multiple genes using the marginal distributions of individual genes, enabling precise capture of data properties and gene expression heterogeneity. In detail, we selected the top 3000 SVGs for each real dataset by SPARK and then generated simulated data by scDesign3 using the filtered gene expression matrix and layer labels as inputs, with the Gaussian process smoother parameter K = 300. Based on the simulated data, we established four scenarios to evaluate the robustness of GraphPCA and other dimension reduction algorithms with varying sequencing depths, noise levels, spot sparsity, and expression dropout rates. Spatial domain detection performance of each method was then evaluated using adjusted Rand index (ARI), normalized mutual information (NMI), and homogeneity score (HS).

GraphPCA model

The GraphPCA model aims to reduce the dimensionality of the gene expression data $X$ while incorporating graph regularization constraints on the low-dimensional embedding. We first recalled the classical PCA, which strives to find a $k$-dimensional projections $Z\in {R}^{n\times k}$, which minimizes the reconstruction error of gene expression matrix $X$. The variational formulation of PCA is given as follows:

$$\underset{Z,W}{min} \, ||X-Z{W}^{\intercal }|{|}_{2}^{2},s.t.{W}^{\intercal }W=I,$$

(1)

where $W$ is an orthonormal eigenvector matrix. This is an optimization problem that can be addressed via the Lagrange multiplier method. The analytical solution to the optimization problem is $Z=XW$, where $W$ represents the top $k$ eigenvectors of ${X}^{\intercal }X$, where $k$ is the dimensionality of $Z$. To ensure the low-dimensional embeddings retain spot location information, a natural idea is to impose graph-based constraints so that physically proximate spots have similar projections. Therefore, we considered the following objective function:

$$\underset{Z,W}{min} \, ||X-Z{W}^{\intercal }|{|}_{2}^{2}+\frac{1}{2}\lambda {\sum }_{\left\{i,j\in A\right\}}{A}_{ij}||{Z}_{i}-{Z}_{j}|{|}_{2}^{2},s.t.{W}^{\intercal }W=I,$$

(2)

where $\lambda$ is a hyperparameter that balances reconstruction error and the smoothness of the projections over the graph. When $\lambda =0$, the GraphPCA degenerates to classical PCA. We noted that the first term along with the constraint, corresponds to the objective of the standard PCA, and the second term is a graph regularization that encourages nearby spot pairs in the graph to have similar projections. The above objective function can be rewritten as

$$\underset{Z,W}{min} \, ||X-Z{W}^{\intercal }|{|}_{2}^{2}+\lambda tr({Z}^{\intercal }LZ),s.t.{W}^{\intercal }W=I,$$

(3)

where $L=D-A$ is the Laplacian matrix of the spatial neighborhood graph, and $D$ is the diagonal matrix of graph $A$. Similar to PCA, the objective function is also non-convex, but we could still derive an analytical solution via the Lagrange multiplier method to solve the $Z$ and $W$ (see Additional file 1: Supplementary Material for details). Let $J$ denote the objective function as given in Eq. (3). We can derive the optimal solution by setting the derivative of $J$ with respect to $Z$ to be zero. Since.

$$\frac{\partial J}{\partial Z}=-2XW+2Z+2\lambda LZ$$

Thus, the derivative equals to zero if and only if ${Z}^{*}={\left(I+\lambda L\right)}^{-1}XW$. To simplify the notation, we denoted

$$K={\left(I+\lambda L\right)}^{-1},$$

which is a symmetric and positive definite matrix. Then, we substituted $Z$ in objective $J$ with ${Z}^{*}=KXW$ and reduced the optimization problem as

$$\underset{W}{min} \, ||X-KXW{W}^{\intercal }|{|}_{2}^{2}+\lambda tr({W}^{\intercal }{X}^{\intercal }{K}^{\intercal }LKXW),s.t.{W}^{\intercal }W=I.$$

(4)

Equation (4) is equivalent to

$$\underset{W}{min} \, tr\left(-{W}^{\intercal }{X}^{\intercal }KXW\right),s.t.{W}^{\intercal }W=I.$$

It is easy to show that the matrix ${X}^{\intercal }KX$ is symmetric and positive semi-definite. The optimal solution ${W}^{*}$ of Eq. (4) is the combination of eigenvectors, associated with the largest $k$ eigenvalues of the graph-revised covariance matrix ${X}^{\intercal }KX$. $X$ corresponding to the top $k$ eigenvalues. It is noteworthy that the solution to GraphPCA can be interpreted as a ridge regularization of the classical PCA, regularized by the Laplacian matrix $L$. Thanks to this closed-form solution, GraphPCA is easy to implement and generates low-dimensional representations in linear time, enhancing integration performance into downstream analyses. As GraphPCA essentially imposes a graph prior on projections, it can be extended to dimensionality reduction tasks for any data with an underlying graph structure.

GraphPCA_multi model

We also developed a multi-sample integration extension called GraphPCA_multi. Suppose a spatial transcriptomics study measures L tissue slices. We followed the same preprocessing procedures as GraphPCA, taking the intersection of the SVGs of all the input slices. For each slice $l=1,\;\cdot\cdot\cdot,\;L$, we constructed a spatial neighborhood graph A_l, with $l=1,\;\cdot\cdot\cdot,\;L$, respectively. Then, we assembled these graphs into a large block diagonal matrix

$${A}_{multi}=diag\left({A}_{1},{A}_{2},\cdots , {A}_{L}\right),$$

where ${A}_{multi}$ represents the integrated spot neighborhood graph. Then, we concatenated the gene expression matrices from all slices to form the input to GraphPCA_multi. Therefore, GraphPCA_multi imposes within-sample smoothness constraints through the block diagonal Laplacian matrix while identifying shared and distinct latent patterns across the multi-slice dataset.

Denoising process

The core principle of the PCA-based framework for denoising lies in the observation that noise typically resides in the lower eigenvalues, while the signal is concentrated in the higher eigenvalues. In the context of GraphPCA, the incorporation of spatial information allows for an enhanced embedding that retains the structural characteristics of the data. This embedding can effectively reconstruct the original gene expression matrix, thereby achieving denoising. In detail, the denoised gene expression matrix $\widetilde{X}$ is computed as the product of the embedding matrix $Z$ and the projection matrix $W$. This approach not only filters out noise but also preserves significant biological signals inherent in ST data.

Downstream analyses

Clustering

We employed the K-means clustering algorithm to detecting spatial domains based on low-dimensional embedding obtained from GraphPCA, PCA, and NMF. Other methods used default embedded clustering algorithms, e.g., SpatialPCA uses a Walktrap algorithm and Louvain algorithm. For fair comparison, the true number of spatial domains was provided to all algorithms. We noted the radius hyperparameter $r$ in DR-SC, defining the neighborhood radius, strongly influences embedding and clustering quality. Therefore, we generated a range of DR-SC embeddings ($r$ from 5 to 100) and chose the best performing for comparisons as the oracle DR-SC version.

Differential gene expression analysis and spatial domain annotation

We employed the Wilcoxon test implemented in the sc.tl.rank_genes_groups function of the Scanpy package to identify differentially expressed genes for each spatial domain with a 1% FDR threshold (Benjamin-Hochberg adjustment). Then, spatial domains are annotated by marker genes and comparing expression spatial patterns against manual annotations.

Trajectory inference

After obtaining the clustering labels, we employed the Slingshot algorithm, the trajectory inference algorithm in scRNA data analysis [82], on the low-dimensional embeddings to depict the spatial trajectory among locations on the tissue, setting white matter and layer 1 as start clusters for DLPFC and mPFC, respectively.

Gene enrichment analysis

We performed gene set enrichment analysis (GSEA [83]) on the top 100 differentially expressed genes sorted by adjusted p values using enrichGO function in the clusterProfiler [84] package, showing the top 30 enriched pathways. Gene sets are downloaded from the Molecular Signatures Database (MSigDB [85, 86], Broad Institute) including C2 (KEGG [87]) and C5 (GO BP: biological process, GO CC: cellular component, GO MF: molecular function).

Clustering performance evaluation metrics

We adopted four different metrics, i.e., adjusted Rand index (ARI), normalized mutual information (NMI), homogeneity score (HS), and McFadden-adjusted pseudo-R² [88], to measure clustering performance of different methods in spatial domain detection. Specifically,

• (1) Adjusted Rand index (ARI)

Given two sets of clustering labels, the ARI is calculated as:

$$ARI= \frac{\sum_{i,j}\left(\genfrac{}{}{0pt}{}{{n}_{ij}}{2}\right)-\left({\sum }_{i}\left(\genfrac{}{}{0pt}{}{{n}_{i}}{2}\right){\sum }_{j}\left(\genfrac{}{}{0pt}{}{{n}_{j}}{2}\right)\right)/ \left(\genfrac{}{}{0pt}{}{n}{2}\right)}{\frac{1}{2}\left({\sum }_{i}\left(\genfrac{}{}{0pt}{}{{n}_{i}}{2}\right)+{\sum }_{j}\left(\genfrac{}{}{0pt}{}{{n}_{j}}{2}\right)\right)-\left({\sum }_{i}\left(\genfrac{}{}{0pt}{}{{n}_{i}}{2}\right){\sum }_{j}\left(\genfrac{}{}{0pt}{}{{n}_{j}}{2}\right)\right)/\left(\genfrac{}{}{0pt}{}{n}{2}\right)},$$

where ${n}_{ij}$ is the number of spots overlapped by cluster $i$ and cluster $j$. ${n}_{i}$ and ${n}_{j}$ are the number of spots in cluster $i$ and $j$, respectively.

• (2) Normalized mutual information (NMI)

NMI is calculated as:

$$NMI= \frac{\sum_{i,j}{p}_{ij}\mathit{log}\frac{{p}_{ij}}{{p}_{i}{p}_{j}}}{\left(-\sum_{i}{p}_{i}\mathit{log}{p}_{i}-\sum_{j}{p}_{j}\mathit{log}{p}_{j}\right)/ 2},$$

where ${p}_{ij}=\frac{{n}_{ij}}{n}$, ${p}_{i}=\frac{{n}_{i}}{n}$, and ${p}_{j}=\frac{{n}_{j}}{n}$.

• (3) Homogeneity score (HS)

The homogeneity score measures the homogeneity or purity within each cluster of the clustering results, calculated as:

$$h=1-\frac{H\left(C|K\right)}{H\left(C\right)},$$

where C represents the ground truth. $K$ represents the cluster labels predicted by the algorithm. $H(C|K)$ is the conditional entropy of the class distribution given cluster labels averaged over clusters, weighted by cluster size:

$$H\left(C|K\right)=-{\sum }_{C=1}^{\left|C\right|}{\sum }_{K=1}^{\left|K\right|}\frac{{n}_{c,k}}{n}log(\frac{{n}_{c,k}}{{n}_{k}}),$$

where ${n}_{c,k}$ is the size of cluster $c$ spots assigned to cluster $k$, ${n}_{k}$ is the size of spots in cluster $k$, and $n$ is the total size of spots. $H(C)$ is the entropy of the class distribution:

$$H\left(C\right)=-\sum_{c=1}^{\left|C\right|}\frac{{n}_{c}}{n}log(\frac{{n}_{c}}{n}).$$

The homogeneity score ranges from 0 to 1, where HS = 1 indicates complete homogeneity, that each cluster only contains spots of a single class.

• (4) McFadden-adjusted pseudo-R² 

The McFadden-adjusted pseudo-R² first constructs a multinomial regression model using low-dimensional features as predictors and the true spatial domains as the outcome, and then calculates the log-likelihood ratios between the regression model and the null model, where no predictors are used. Within the regression model, we calculated the McFadden-adjusted pseudo-R² to assess the predictive capacity of the predictor variables in projecting the ground truth. A higher pseudo-R² indicates that the method can effectively extract informative output for predicting the true spatial domains.