# jsPCA: fast, scalable, and interpretable identification of spatial domains and variable genes across multi-slice and multi-sample spatial transcriptomics data

Ines Assali, Paul Escande, Paul Villoutreix

► To cite this version:

Ines Assali, Paul Escande, Paul Villoutreix. jsPCA: fast, scalable, and interpretable identification of spatial domains and variable genes across multi-slice and multi-sample spatial transcriptomics data. 2025. hal-05391223

HAL Id: hal-05391223

https://hal.science/hal-05391223v1

Preprint submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# jsPCA: fast, scalable, and interpretable identification of spatial domains and variable genes across multi-slice and multi-sample spatial transcriptomics data

Ines ASSALI<fnref n="1" />, Paul ESCANDE<fnref n="2" /> and Paul VILLOUTREIX<fnref n="1" />

Corresponding author: ines.assali@univ-amu.fr, paul.villoutreix@univ-amu.fr

**Keywords** Spatial transcriptomics, Spatial statistics, Machine learning, Spatial Domains, Spatially Variable Genes

**Abstract** Spatially structured cell heterogeneity within tissues is essential for healthy organ function. This heterogeneity is reflected by differential gene expression activity at various spatial location. Spatial transcriptomics technologies record genome-wide measurements of gene expression at the scale of entire tissues with high spatial resolution. While they have revolutionized our quantitative understanding of tissue architecture, these technologies generate large and high dimensional datasets encompassing tens of thousands of genes recorded at tens of thousands of spatial locations, requiring efficient automated methods for their analysis. In this study we introduce joint spatial PCA (jsPCA), a novel, fast, scalable and interpretable method for the automatic identification of spatial domains and variable genes in multi-slice and multi-sample spatial transcriptomics data. jsPCA relies on a simple mathematical formulation of a spatial covariance defined as the product of the gene expression covariance with the spatial autocorrelation. The principal components of this spatial covariance yield a biologically meaningful low-dimensional representation. From this representation, we can derive spatial domains by simple clustering. In addition, spatially variable genes can be identified directly from the principal components coefficients. Moreover, this approach enables the joint representation of multiple slices and samples, a frequent experimental setting. This joint representation is obtained without spatial alignment by computing common principal components via joint diagonalization of the set of spatial covariance matrices obtained for each slice. By leveraging sparsity and non-convex optimization on manifold, jsPCA leads to computing time in the order of seconds to minutes, substantially outperforming state-of-the-art approaches. We benchmarked jsPCA on the Visium 10x dataset of human dorsolateral prefrontal cortex and the Stereo-seq MOSTA dataset of mouse embryonic development against 10 state-of-the-art methods. Our approach demonstrated excellent performances, comparable or better than state-of-the-art methods, such as SpatialPCA, BASS, GraphPCA or Stagate, while being much faster, interpretable, and scalable to very large datasets.

## Introduction

Spatial transcriptomics (ST) has changed the way we understand biological systems by offering a perspective that combines gene expression and spatial information within tissues. Unlike traditional single-cell sequencing technologies that dissociate cells from their microenvironment, ST preserves the spatial structure of tissues while capturing gene expression profiles, enabling an in-depth analysis of cellular expression *in situ* [1]. This technological advance has enabled researchers today not only

<footnote n="1">¹ Aix-Marseille Université, MMG, Inserm U1251, Turing Centre for Living systems, Marseille, France</footnote>

<footnote n="2">² Institut de Mathématiques de Toulouse; UMR 5219, Université de Toulouse, CNRS; UPS, F-31062 Toulouse Cedex 9, France</footnote>
to map the location of gene expression in tissues but also to understand the mechanism of cell-cell
interaction [2], to identify cell types according to their localization and molecular profile [3], and to
better understand a variety of biological processes, including tissue pathology [4], immune responses
[5], and embryonic development [6]. These advances open up new possibilities for personalized
medicine and the development of targeted treatments [7]. However, this technology faces significant
challenges in terms of data sparsity, very high dimensionality, and technical and biological noise [8].
These fundamental data limitations require innovative computational strategies capable of efficiently
normalizing, reducing dimensionality while preserving essential biological insights.

Specifically, the task of identifying spatial domains that are both consistent in terms of expression and
spatial organization has gathered attention in the past few years [9]. Several approaches have been
proposed to solve this task as a problem of spatially aware clustering. Some of these methods rely on
graph neural networks, such as STAGATE, which is a graph attention auto-encoder framework [10],
SpaceFlow, which is a deep learning framework based on spatially regularized deep graph networks
[11], and CCST (Cell Clustering for Spatial Transcriptomics), which is an unsupervised cell clustering
method based on a graphical neural network [12]. Another set of methods relies on a Bayesian frame-
work, such as BASS which is a multi-scale, multi-sample approach to spatial transcriptomic analysis
[13]. BASS explores data at different scales by simultaneously clustering cell types and detecting
spatial domains. These analyses are integrated into a hierarchical Bayesian modeling framework.
Similarly, BayesSpace is a Bayesian statistical method that uses the neighborhood information [14].
Finally, some approaches have been proposed to extend the classical Principal Component Analysis
approach to spatial transcriptomics, such as SpatialPCA [15], which is based on probabilistic PCA, it
includes location information as an additional input and uses a kernel matrix to model spatial corre-
lations between tissue sites. Similarly, GraphPCA is a dimension reduction algorithm that integrates
a graph based regularization of principal component analysis (PCA) to improve the analysis of spatial
transcriptomic data [16]. These methods are complex, not always interpretable, or require tuning
several hyperparameters, and can be slow in their execution.

In addition to spatial domain identification, an important step in spatial transcriptomics analysis is
identifying which genes are spatially variable (SVGs) [17]. These genes exhibit non-random expres-
sion patterns across tissue, often reflecting specific biological functions or distinct cellular microen-
vironments. Their identification not only highlights relevant spatial patterns, but can also be used
to reduce the dimensionality of the data. SVGs are similar to the selection of highly variable genes
(HVGs) in single-cell transcriptomics; however, they differ in that they preserve the spatial relation-
ships specific to biological samples. SVGs thus play a central role in many downstream analyses,
such as the exploration of pathological mechanisms [17]. Numerous approaches have been devel-
oped for this task, ranging from methods based on Gaussian processes (SpatialDE) [18] to mixed or
non-parametric models (SPARK [19], SPARK-X [20]), including self-organizing maps (SOMDE) [21], sta-
tistical enrichment in spatial networks (Giotto) [22], Gaussian processes based on nearest neighbors
(NnSVG) [23], and spatial autocorrelation measures such as Moran's I (implemented in Seurat [24],
among others These approaches have mostly been developed independently of domain identifica-
tion methods and are usually used as a preprocessing step to reduce the dimensionality of the data.
Moreover, they rely on sometimes complicated statistical models.
In this study, we propose joint spatial Principal Component Analysis (jsPCA), a novel computational method for identifying spatial domains and spatially variable genes from spatial transcriptomics data that addresses the limitations of existing methods. jsPCA generalizes sPCA, an established method for population genetics previously applied to genetic variability using allelic frequency data and spatial coordinates [25]. We previously introduced an initial formulation of jsPCA to study single-cell morphometrics data within biological tissues obtained by microscopy imaging [26]. Adapting jsPCA to spatial transcriptomics data generates new methodological challenges related to the high dimensionality of this type of data. We solved these by exploiting the sparse nature of the matrices and non-convex optimization on manifold techniques [27,28]. jsPCA efficiently integrates spatial information with gene expression profiles to produce a low-dimensional representation of the data that enable the identification of spatial domains and spatially variable genes. jsPCA is well adapted for multi-slice and multi-sample data.

# Methods

## Data preprocessing

To optimize the detection of biological signals, we carried out a three steps pre-processing of the gene expression data using the Scanpy library [29]. Firstly, we applied a prevalence-based selection criterion, selecting only genes expressed in at least 20 cells/spots. This eliminates potentially artefactual and poorly informative variables. Next, we transformed the data using Pearson's residual normalization, which is a technique that scales expression levels based on the relationship between mean and variance, allowing us to correct any technical biases while preserving biologically significant variability. Then we normalized the data, transforming each gene expression profile into z-scored values, guaranteeing zero mean and unit variance. This transformation makes expression levels across different genes comparable, regardless of their initial expression amplitude, facilitating the robust identification of relevant biological patterns. Normalization is performed separately for each slice and sample.

## Construction of spatial neighborhood graph

To encode the spatial information contained in ST, we built a neighborhood network based on the spatial coordinates of the spots/cells. For this, we applied the *k*-nearest neighbor (*k*-NN) algorithm based on Euclidean distance, adapting the classic elbow method [30] to determine the optimal number of *k* neighbors. This modified version selects *k* by analyzing the variation in average distances to the *k*-th neighbors and identifying the point where the slope of the distance curve changes most strongly. The connectivity of the resulting network is encoded into a sparse adjacency matrix *C*. The rows and columns of *C* correspond to the spots/cells and $C(i, j) = 1$ if $j \in N(i)$, $C(i, j) = 0$ otherwise. *N(i)* represents the set of neighbors of spot *i*. We denote by *L* the row-normalized version of the matrix *C*, such that the sum of each row is equal to one.

## Join spatial principal component analysis (jsPCA)

jsPCA combines principal component analysis and spatial autocorrelation for multivariate spatial variables and multiple samples. Fig. 1 shows a global overview of the method. It is based on the idea of finding an embedding of the data that represents best both the covariance, as in a classical PCA approach, and the spatial autocorrelation (also known as Moran's index) derived from the
spatial information [25,31,26]. In the context of spatial transcriptomics data, the problem consists in identifying gene expression patterns while preserving the spatial organization of tissues. Transcriptomic covariance is calculated using the $n \times p$ gene expression matrix $X$, and the Moran's index is calculated as $X^T LX(X^T X)^{-1}$.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Overview of the Joint spatial PCA (jsPCA) approach. jsPCA is a joint spatial principal component analysis method that simultaneously integrates variation in gene expression and spatial structure in spatial transcriptomics (ST) data. It starts with a ST dataset consisting of a gene expression matrix $X_j$ and the spatial coordinates of tissue cells/spots. From these coordinates, we construct a spatial neighborhood graph, represented by a normalized adjacency matrix $L_j$, which encodes the relationships between spots. For each dataset, jsPCA then calculates a spatial covariance matrix $A_j$ by combining gene expression covariance and spatial autocorrelation (measured by Moran's index). A joint eigenvalue analysis is then performed on the set $\{A_1, \dots, A_m\}$, producing a shared basis of eigenvectors $P$ describing the main spatio-molecular patterns. Each dataset is projected into this latent space via $X_j P$, and these representations are then used to group spots into coherent spatial domains within the tissue while the principal components can be directly interpreted in terms of spatially variable genes.</figcaption></figure>

For a given sample, with expression data $X$ and normalized adjacency matrix $L$, representing best both the covariance and the spatial autocorrelation can be modeled as the question of identifying the main directions of the product matrix between these two quantities: $X^T LX(X^T X)^{-1} \times \frac{1}{n} X^T X = \frac{1}{n} X^T LX$. It was previously shown in [25] that the principal components (PCs) of this product matrix can be found as the eigenvectors $P$ with eigenvalues $D$ of the symmetrized matrix, which we call spatial covariance matrix:

$$A = \frac{1}{2n} X^T (L + L^T) X \in \mathbb{R}^{p \times p}. \qquad (1)$$

To tackle multiple datasets, such as multislice or multisample ST datasets, we need to be able to simultaneously process multiple datasets with full consideration of the spatial structure of each individual slice or sample [26]. Given a collection of samples $\{X_j, L_j\}$, jsPCA computes a joint diagonalization of the corresponding spatial covariance matrices, this can be summarized by finding the common set of eigenvectors $P$ such that

$$A_j := \frac{1}{2n_j} X_j^T (L_j + L_j^T) X_j \approx P D_j P^T, \qquad (2)$$
where $D_j$ is the diagonal matrix corresponding to the eigenvalues specific to sample j. We show in the supplementary material how we obtain numerically the joint eigenvectors by non-convex optimization on a manifold because of the large size of the matrices in ST (Supplementary text 1).

The new representation $\tilde{X}_j$ of the initial data $X_j$ is obtained by projecting the data on the eigenvectors $P$. In practice we compute $\tilde{X}_j = X_j P_k$ where $P_k$ represents the top $k$ eigenvectors of the common set of eigenvectors $P$. In the case of the analysis of a unique ST dataset, $\tilde{X}_j$ is obtained equivalently by projecting $X_j$ on the eigenvectors of the spatial covariance $A_j$. In the experiments below, the actual value of $k$ is chosen by cross-validation in a range of 10 to 50.

Several approaches for the computation of spatial domains in spatial transcriptomics have been derived from a formulation of Principal Component Analysis with a spatial component. In the supplementary material (Supplementary text 2), we provide a formal comparison with these four methods (SpatialPCA [15], GraphPCA [16], RASP [32], nichePCA [33]). We show that jsPCA differs by the absence of hyper-parameters controlling the influence of spatial information and the absence of kernel. In contrast to other methods, jsPCA has a very simple and general formulation.

## Spatial Domain Identification

For spatial domain identification, we use the low dimensional representation $\tilde{X}_j$ defined above, on which we apply a simple clustering algorithm, the Gaussian Mixture Model (GMM). In our experiment, the number of domains is given by the ground truth annotation, but could in principle be obtained by minimizing the Akaike Information Criterion or the Bayesian Information Criterion [34]. The labels obtained can then be reported on the tissue and compared to the ground truth.

Following clustering, we applied a refinement step, as is common in most spatial domain identification methods [9]. This refinement step was initially proposed in SPAGCN [35] and consists in a very simple step. For each spatial location, we compare its predicted label to the ones of its direct neighbors. If more than half of them have a different label, then we relabel the current one with the most common one within its direct surrounding.

## Spatially Variable Gene

Our approach naturally leads to a list of spatially variable genes [36]. These genes can be directly obtained as those with the highest coefficients of the top principal component. This strategy enables the rapid and robust identification of SVGs in large-scale datasets. Other approaches such as SPARK-X [20] rely on a statistical method that detects SVGs using a robust covariance test simultaneously integrating gene expression and the spatial organization of cells/spots. Specifically SPARK-X constructs a test statistic based on the trace of the product between the expression and distance covariance matrices and evaluates it with a $\chi^2$ mixture distribution to obtain calibrated p-values for selecting significant genes. jsPCA offers a simpler and more computationally efficient procedure.

## Datasets and Method Evaluation

Our approach was first evaluated using the Visium 10x dataset of human dorsolateral prefrontal cortex (DLPFC) [37]. This dataset was obtained as part of the spatialLIBD project (http://spatial.libd.org/spatialLIBD). The dataset is composed of 12 slices acquired from three different individuals (four slices per sample). The number of genes per slice is 33538, with a median of 3844 spots.
The samples were annotated manually, each sample comprising five or seven layers, including six cortical layers (L1 to L6) and white matter (WM), see Fig. 3a.

In the first scenario, we evaluated the accuracy of spatial domain identification by jsPCA. To do so, we benchmarked jsPCA against STAGATE [10], SpatialPCA [15], GraphPCA [16], BASS [13], CCST [12], BayesSpace [14], Spaceflow [11], Leiden and Louvain [38] on the same dataset. Each spatial clustering method applies a specific pre-processing technique. In this comparison, we have adopted the optimal strategies recommended by the authors of each method. To have a fair comparison with the other methods, which do not handle multi-slice and multi-sample datasets, we performed a monoslice analysis using jsPCA, where we considered each slice of the 12 slices separately. For each slice, we projected the expression matrix on the first 10 to 50 (depending on the slice) eigenvectors then performed clustering on the results obtained.

In the second scenario, we performed a multi-slice analysis with jsPCA, exploring three configurations: 2 slices, 4 slices and 12 slices. In each of these cases, we considered jointly groups of slices, for which we computed a joint embedding embedding with jsPCA and performed clustering simultaneously. In the 2-slice analysis, we examined each pair of nearby slices, separated by a gap of 10 micrometers. In the 4-slice analysis, we studied groups of 4 slices belonging to the same sample with gaps of 10 micrometers, 300 micrometers and 10 micrometers between them. Finally, in the 12-slice analysis, we studied the three samples together, each sample consisting of 4 slices with the same gaps (10 micrometers, 300 micrometers and 10 micrometers).

We also evaluated our approach on the MOSTA (Mouse Organogenesis Spatiotemporal Transcriptomic Atlas) [39] dataset, which consists of 53 sagittal sections of mouse embryos collected at developmental stages from E9.5 to E16.5 and profiled using Stereo-seq technology. The slices contain a median of 84811 spots and 27810 detected genes, providing a detailed representation of transcriptomic organization during development. For a fair comparison, we first carried out a monoslice analysis with jsPCA, in the same manner as with the DLPFC dataset, treating each sagittal section independently, and benchmarked the results against established methods such as STAGATE, SpatialPCA, GraphPCA, and BASS, using the preprocessing strategies recommended by their authors. We then extended the analysis to a multi-slice setting by jointly considering all sections from the same developmental stage, covering E9.5 to E13.5 (5 slices at E9.5, 4 at E10.5, 4 at E11.5, 6 at E12.5, and 4 at E13.5).

Clustering performance was estimated using Adjusted Rand Index (ARI) and Normalized Mutual Information (NMI). To evaluate scalability, we measured the elapsed time and peak memory usage of our approach as well as the other approaches, on a machine with 32 GB of RAM and 24 cores (Ubuntu v.24.04), while all computations for the MOSTA dataset were executed on the IFB (Institut Français de Bioinformatique) national cluster, which provides high-performance computing resources.

## Results

To measure the ability of jsPCA at identifying spatial domains, we first evaluated it and compared it to state-of-the-art methods on the DLPFC reference dataset. We showcase the ground truth (the six manually annotated layers and white matter (WM) of the DLPFC dataset) as well as the domains identified (cluster assignments) by the different methods for sample 151673 (Fig. 2a). This compari-
son demonstrates the well-balanced performance of jsPCA for the identification of spatial domains
in spatial transcriptomics. With an ARI of 0.62 and an NMI of 0.71, jsPCA achieves performances
comparable to, or better than, the best existing approaches, such as SpatialPCA (ARI 0.57, NMI 0.70),
BASS (ARI 0.59, NMI 0.70), and GraphPCA (ARI 0.47, NMI 0.65). Importantly, jsPCA is far more effi-
cient, requiring only 22.7 seconds and 623.7 MB of memory, compared to 235.2 s / 7571.7 MB for
SpatialPCA, 306.1 s / 3504.2 MB for BASS, and 624.4 s / 457.7 MB for GraphPCA. Unlike Bayespace,
with its extremely high memory consumption (8699.9 MB) and longer runtime (893.7 s), or CCST,
with its lower ARI and NMI scores (0.31 and 0.46) and higher memory usage (2137.2 MB), jsPCA of-
fers a good compromise between accuracy and computational efficiency. In addition, our method
outperforms Spaceflow (ARI 0.17, NMI 0.45), which achieves lower performance and produces nois-
ier segmentation. As expected, algorithms that do not take spatial information into account, such as
Leiden (ARI 0.21, NMI 0.31) and Louvain (ARI 0.22, NMI 0.30), fail to distinguish the different cortical
layers, leading to low clustering accuracies. These results underline the importance of effectively
integrating spatial information into dimension reduction methods to optimize domain identification
while preserving a moderate computational cost.

In analyzing each of the sections contained in the entire database (Fig. 2b), we noted that of all the
spatial analysis methods compared in this study, jsPCA stands out for its outstanding balance be-
tween precision and reproducibility (remarkable stability and low variability of results). Its clustering
performance, measured by ARI and NMI, reaches substantial median values, placing it in the top tier
with SpatialPCA and Bass, but well above the other approaches. The real strength of jsPCA lies in its
consistency, illustrated by the narrowness of its interquartile range for ARI and NMI. This low disper-
sion indicates that results remain reliable across different slices, unlike methods such as SpatialPCA
(which has outliers) or STAGATE (which has higher variability). The total absence of outliers for jsPCA
confirms its robustness in the face of special cases that could degrade the performance of other
methods. Comparative analysis (Fig. 2c) reveal that jsPCA also stands out for its computational effi-
ciency, combining one of the shortest execution times (around 50 seconds) and moderate memory
usage (800 MB), far below the requirements of SpatialPCA (5000 MB) or Bayespace (8000 MB).

After finding that, in its simplest setting, i.e. jsPCA on a single ST slice, our approach performs in a
very fast and efficient way at accuracy level better or comparable to the top methods of the state
of the art, we inquired the impact of considering multiple slices at the same time with jsPCA. Con-
sidering multiple slices at the same time has the advantage of predicting domains that are directly
comparable between slices, as they are computed from a joint representation. Fig. 3b shows the
evolution of clustering performance, evaluated using the NMI and ARI indices, as a function of the
number of slices. Accuracy is analyzed for different configurations: monoslice, 2-slice, 4-slice and 12-
slice. While monoslice remains the most accurate, the 2-slice and 4-slice configurations still achieve
good performance with only modest reductions. In contrast, the 12-slice configuration shows lower
accuracy, which is expected since it merges slices coming from three different DLPFC samples, and
these cross-sample differences reduce the consistency of the clustering. In terms of computational
performances, however, there is a systematic increase of time and memory when increasing the
number of slices considered together in jsPCA (Fig. 3b).
<figure><img src="image_3.png" /><figcaption>**Fig. 2. Performance of jsPCA for spatial domain identification in the human dorsolateral prefrontal cortex (DLPFC) 10x Visium dataset.** (a) Spatial domain prediction results produced by different methods for slice 151673 of the DLPFC dataset, shown alongside the ground truth annotation of the six cortical layers and white matter (WM). ARI, NMI, runtime, and memory usage values are reported for each method. (b) Box plots showing clustering accuracy for all 12 DLPFC sections, measured by NMI (left) and ARI (right) across all methods. (c) Comparison of runtimes (left) and memory usage (right) across all methods.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** jsPCA performance across different configurations (monoslice, 2-slice, 4-slice and 12-slice) on the DLPFC dataset. (a) Distribution of the slices among three samples in the DLPFC dataset (b) Comparison of NMI and ARI scores and evaluation of runtime and memory usage.</figcaption></figure>

jsPCA yields a representation of the data, the eigenvectors, that can be used independently of the dataset on which they have been computed. Therefore, jsPCA can in principle be used to predict a good embedding of an unseen dataset. To evaluate the prediction power of jsPCA, we performed a cross-slice analysis on the DLPFC dataset (Fig. 4a). For each sample, we trained jsPCA on three slices, learning the joint eigenvectors (principal components) and clustering parameters from the GMM model, and then used them to predict the domain labels of the remaining held-out slice. This process was repeated in a leave-one-out manner so that every slice served once as a test case. This approach makes it possible to test how well jsPCA generalizes from one slice to another within the same sample, showing whether it can handle structural differences while still producing consistent and reliable predictions. As shown in (Fig. 4b), the monoslice configuration reaches somewhat higher NMI score, but the cross-slice setting delivers strong and competitive results and even improve the ARI score. This indicates that jsPCA is able to generalize effectively across slices, handling structural variability without a major loss in accuracy.

In Fig. 5a, we compare the performance of jsPCA with several reference methods (STAGATE, SpatialPCA, GraphPCA, and BASS) for identifying spatial domains in the MOSTA dataset (E11.5_ES1). The domains obtained by jsPCA show a high degree of agreement with the manual annotation and achieves the best quantitative scores, with an NMI of 0.67 and an ARI of 0.48, outperforming the other methods. In Fig. 5b, we present a monoslice analysis covering slices from developmental stage E9.5 through the first slice of stage E12.5. This subset was chosen because SpatialPCA and GraphPCA could not process the remaining slices of the dataset due to excessive memory requirements. Quantitative comparisons confirm the robustness of jsPCA, which achieves a strong balance between clustering accuracy and stability. The evaluation of runtime and memory usage demonstrates that jsPCA is computationally efficient relative to alternative methods, requiring on average only around 97 seconds and around 8.1 GB of memory, compared to several thousand seconds and much higher
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Cross-slice evaluation of jsPCA on the DLPFC dataset. (a) Schematic illustration of the cross-slice analysis. Each sample contains four slices. jsPCA was trained on three slices (blue boxes) to learn joint embeddings and clustering, and then used to predict the domains of the held-out slice (brown box). This procedure was repeated in a leave-one-out manner so that each slice served once as test data. (b) Quantitative comparison of monoslice (blue) and cross-slice (brown) configurations, measured by Normalized Mutual Information (NMI) and Adjusted Rand Index (ARI).</figcaption></figure>
memory usage for SpatialPCA, GraphPCA, and BASS. In Fig. 5c, we extend the monoslice analysis to include all slices of the dataset. jsPCA maintains strong clustering performance in terms of NMI and ARI, with stable results across the full range of slices. When considering runtime, jsPCA continues to be highly efficient, completing analyses much faster than STAGATE or BASS. In terms of memory, the requirements increase compared to the subset analysis, which is expected since the remaining slices contain a substantially larger number of spots. This higher memory demand explains why SpatialPCA and GraphPCA could not process these slices, whereas jsPCA successfully scaled to the entire dataset while preserving accuracy and robustness.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Performance of jsPCA for spatial domain identification in the mouse organogenesis spatio-temporal transcriptomic atlas (MOSTA) Stereo-seq dataset. (a) Ground truth obtained from [39]. Spatial domain identification results from the jsPCA, STAGATE, SpatialPCA, GraphPCA, and BASS methods on the E11.5.E1S1 sagittal section, shown alongside the manual reference annotation of major organs and tissues. ARI, NMI, runtime, and memory usage are reported for each method. (b) Comparative analysis on the first 14 sagittal sections of MOSTA, with box plots showing NMI (left) and ARI (right) scores for jsPCA, BASS, STAGATE, GraphPCA, and SpatialPCA, along with their runtimes and memory usage. (c) Additional evaluation on all sagittal sections of MOSTA for jsPCA, BASS, and STAGATE, using the same accuracy, runtime, and memory usage metrics.</figcaption></figure>
Fig. 6 compares the performance of jsPCA on the MOSTA dataset in monoslice and multi-slice configurations. As shown in Fig. 6b, the monoslice configuration has overall slightly higher NMI and ARI scores, while the multi-slice configuration remains competitive, showing how robust jsPCA is in both cases. It also shows that analyzing multiple slices together requires more time and memory, which is predictable, but this increase in computational cost does not come at the expense of clustering quality. These results highlight jsPCA's ability to deliver reliable and stable performance across different analysis settings.

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** MOSTA performance across different configurations (monoslice, multislice). (a) jsPCA monoslice spatial domain identification across all the time points. (b) Comparison of NMI and ARI scores and evaluation of runtime and memory usage between monoslice and multislice settings on 5 slices at E9.5, 4 at E10.5, 4 at E11.5, 6 at E12.5, and 4 at E13.5.</figcaption></figure>

Finally, to evaluate the ability of jsPCA to identify spatially variable genes (SVGs), we compared its outputs with those obtained from SPARK-X. The results show that jsPCA captures a substantial proportion of the SVGs identified by SPARK-X, particularly in its first component (denoted jsPCA1). In the DLPFC dataset (Fig 7a), jsPCA1 recovers about 82% of the SVGs detected by SPARK-X, while in the MOSTA dataset (Fig 7b) the overlap is nearly 79%. The distribution of jsPCA coefficients shows that the top-ranked genes contribute most of the spatial signal, which justifies selecting the top 3000 SVGs. When considering the top 3000 SVGs ranked by jsPCA1, the vast majority of SPARK-X genes are consistently present. These results highlight that jsPCA not only provides a simpler and more computationally efficient alternative to SPARK-X but also achieves a high level of agreement in the genes it identifies, supporting its reliability for large-scale spatial transcriptomics analyses.

## Discussion

We introduced jsPCA, a computational framework for efficient and interpretable analysis of spatial transcriptomics data. By integrating gene expression covariance with spatial autocorrelation, jsPCA identifies both spatial domains and variable genes in a single parameter-free formulation. Com-
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Detection of spatially variable genes by jsPCA: comparison with SPARK-X on (a) DLPFC (slice 151673) and (b) MOSTA (E16.15_E1S1). In the first column, we report the percentage of genes common between the 3000 most spatially variable genes identified by SPARK-X and the 3000 genes with the highest coefficients in each of top 10 eigenvectors. In the second column, we show the distribution of the coefficient in the top jsPCA eigenvector, the red dashed line represents the 3000th top coefficient. In the third column, we show how the genes common between SPARK-X and jsPCA are distributed in the jsPCA eigenvector, the x-axis represents the rank of the gene according to its coefficient value, the plots indicate if the gene is present or not in the SVGs identified by SPARK-X.</figcaption></figure>

pared with existing approaches, jsPCA consistently provides a favorable trade-off between accuracy, interpretability, and computational efficiency.

A key contribution of jsPCA is its ability to jointly analyze multiple tissue slices or samples without explicit spatial alignment. This contrasts with existing PCA-inspired methods, which typically operate on single slices or rely on distance kernels that do not generalize well across samples. We show that jsPCA generalizes across slices within the same tissue by performing accurate domain prediction on unseen slices. Moreover, jsPCA can be applied jointly to multiple samples, opening the way to integrative analyses of multi-sample ST datasets.

Our results also demonstrate that jsPCA is scalable to atlas-level datasets. On the Stereo-seq MOSTA dataset, jsPCA handled dozens of slices with tens of thousands of spots, whereas several leading methods failed due to excessive memory requirements. With the rapid emergence of ultra-high-resolution ST platforms such as Visium HD and whole-organ Stereo-seq, scalability is becoming a practical bottleneck; jsPCA offers a solution by maintaining efficiency without compromising interpretability.

Beyond benchmarking, jsPCA provides a framework for biological discovery. Because the principal components directly highlight spatially variable genes, jsPCA yields interpretable axes of variation that can be linked to known or novel tissue structures. For example, in the mouse embryonic dataset, jsPCA not only reproduced manual annotations of organ primordia but also identified subtle subdomains that were not detected by competing methods. This underscores the potential of jsPCA to uncover new biology in atlas-scale datasets.
Nevertheless, some limitations remain. jsPCA does not explicitly correct for batch effects, which can confound multi-sample comparisons. However, the modular nature of the framework allows for integration with established batch correction tools upstream of analysis [40]. Another limitation is that while jsPCA naturally extends to any spot- or cell-based dataset [26], its applicability to truly multi-modal omics data remains to be established. Future work will aim to extend jsPCA to joint spatial analysis across transcriptomic, proteomic, and metabolomic layers [41].

Overall, jsPCA provides an accessible, scalable, and biologically interpretable framework for spatial domain and gene identification in spatial transcriptomics. By enabling multi-sample joint analysis and scaling to datasets that currently exceed the reach of competing methods, jsPCA is positioned to become a standard tool for large-scale spatial omics studies, particularly in developmental biology and tissue atlas projects.

## Availability and Implementation

An open-source implementation of the jsPCA algorithm can be found on https://github.com/VILLOUTREIXLab/jsPCA.

## Acknowledgments

We thank Solène Song, Denis Puthier, Anaïs Baudot, Franck Picard for fruitful discussions. This research was supported by Agence Nationale de la Recherche (ANR) project ANR-CPJ12022-VILLOUTREIX-U125. It benefited from access to the Core Cluster of the Institut Français de Bioinformatique (IFB).

## Supplementary Text 1 - Computation of jsPCA for large matrices

To compute jsPCA eigenvectors for a collection of ST datasets, i.e. a collection of expression matrices and normalized adjacency matrices {$X_j, L_j$}, we need to perform a joint diagonalization of the corresponding spatial covariance matrices. This can be summarized by finding the common set of eigenvectors $P$ such that

$$A_j := \frac{1}{2n_j} X_j^T (L_j + L_j^T) X_j \approx P D_j P^T,$$

The matrices $P$ and {$D_j$} can be obtained by solving

$$ \begin{array}{ll} \min_{P,D_j} & \frac{1}{2} \sum_j \|PD_jP^T - A_j\|^2 \\ \text{subject to} & P \in \mathbb{R}^{p \times p}, \quad P^T P = PP^T = I_p \\ & \lambda_j \in \mathbb{R}^p, \quad D_j = \text{diag}(\lambda_j) \end{array} $$

Given the large size of the matrices handled in ST, the joint diagonalization of the spatial covariances matrices {$A_j$} cannot be obtained by classical approaches such as Jacobi-like techniques because they require too much memory [42]. Instead, we implemented the following approach.
## Independent truncated diagonalizations

First, the truncated diagonalizations of order $K$ of each $A_j$ are computed as $A_j \approx P_j S_j P_j^T$. In our experiments, we use values of $K$ ranging from 10 to 50, hence $K \ll p$. For $K \ll p$ the truncated diagonalization can be efficiently computed using an Implicitly Restarted Lanczos Method which is implemented in the scipy library [43]. The key feature of this method is its *matrix-free* operation meaning that it only uses matrix-vector products with $A_j$ and does not require a matrix built in memory. However, since $p$ is large and the matrices $A_j$ are dense, directly computing the eigensystem from $A_j$ is particularly inefficient (considering the building time and the dense matrix-vector products). Fortunately the matrices $X_j$ and $L_j$ are very sparse. Implementing an efficient matrix-vector product with $A_j$ can therefore be performed by successive sparse matrix-vector products, that is implementing $A_j u = \frac{1}{2n_j} X_j((L_j + L_j^T)(X_j^T u))$ for $u \in \mathbb{R}^p$.

## Joint approximate diagonalization

Given the eigensystems $\{(P_j, D_j)\}_{j}$, the joint diagonalization can now be modeled as

$$ \begin{array}{ll} \min_{P,D_j} & \frac{1}{2} \sum_j \|PD_jP^T - P_jS_jP_j^T\|^2 \\ \text{subject to} & P \in \mathbb{R}^{p \times K}, \quad P^T P = I_K \\ & \lambda_j \in \mathbb{R}^K, \quad D_j = \text{diag}(\lambda_j) \in \mathbb{R}^{K \times K} \end{array} $$

With an argument similar to [44] this problem can be shown to be equivalent to

$$ \begin{array}{ll} \min_{P} & -\frac{1}{2} \sum_{j} \| \text{diag}(P^T P_j D_j P_j^T P) \|^2 \\ \text{subject to} & P \in \mathbb{R}^{p \times K}, \quad P^T P = I_K \end{array} $$

which can be solved by a gradient descent on the Stiefel manifold [27,28] (modeling the orthogonality constraint $P^T P = I_K$). Following again [44], an initialization close to a minimizer can be obtained by computing the $K$ largest left singular vectors of the matrix $(P_1 S_1, \dots, P_j S_j, \dots)$ (i.e. concatenating the columns of the $\{P_j S_j\}$) and by letting $D_j = \text{diag}(P^T P_j S_j P_j^T P)$. In our experiments, we performed 100 steps for the gradient descent.

## Supplementary Text 2 - Formal comparison between jsPCA and spatial-PCA like approaches

We review in this section the differences between various formulations of spatial PCA applied to spatial transcriptomics that have been proposed recently, with harmonized notations.

We consider a matrix of gene expression measurements for $p$ genes in $n$ spatial locations of a tissue, denoted by $X \in \mathbb{R}^{n \times p}$. These measurements have spatial coordinates that are denoted by $S = \{s_i \in \mathbb{R}^k, i \in \{1, \dots, n\}\}$ where $k = 2$ or 3.

Recall the variational formulation of PCA, which consists in finding a $d$-dimensional projection $Z \in \mathbb{R}^{n \times d}$ such that:
$$ \min_{Z,W} \|X - ZW^T\|_2^2, \text{ with } W^T W = I \qquad (3) $$

This optimization problem can be solved analytically. The analytical solution leads to $Z = XW$ where $W$ represents the top-$d$ eigenvectors of the empirical covariance matrix $X^T X$.

PCA however does not take into account the spatial information $S$. Several methods have been proposed to incorporate spatial information in the context of spatial transcriptomics data analysis. We will distinguish between two families of approaches, the approaches based on the distance between spots and the approaches based on a spatial graph. In that case, the locations are turned into an undirected graph with the $k$-nearest neighbor algorithm, we denote by $C$ the corresponding adjacency matrix. $C_{i,j} = 1$ if $j$ is in the neighborhood of $i$, 0 otherwise.

## Distance-based spatial PCA

**Spatially aware dimension reduction for spatial transcriptomics [15]** The expression matrix $X$ is assumed to have been normalized through variance stabilizing transformation and further scaled for each gene to have zero mean and unit standard deviation.

The goal is to infer a $n \times d$ factor matrix $Z$ that represents a low dimensional embedding of $X$. The approach is based on the following latent factor model:

$$ X = C_v B + ZW^T + E \qquad (4) $$

where $C_v$ is an $n \times q$ matrix of covariates, and $B$ is a $q \times p$ matrix of corresponding coefficients, here $q$ is assumed to be equal to 1. $W$ is a $p \times d$ factor loading matrix and $E$ is an $n \times p$ matrix of residual errors, with elements $E_{i,j}$ following independent normal distribution with mean zero and variance $\sigma_0^2$, $E_{i,j} \sim N(0, \sigma_0^2)$.

$W$ is assumed to have orthonormal columns such that $W^T W = I$. For the $d$ factors $Z_l$, i.e. the vectors of dimension $n$ composing the factor $Z$ matrix, the assumption is that they follow a multivariate normal distribution:

$$ Z_l \sim N(0, \Sigma_l) \qquad (5) $$

with $\Sigma_l$, the $n \times n$ covariance matrix modeling the correlation among the spatial locations. $\Sigma_l$ has for elements $\Sigma_l(i, i') = \sigma_0^2 \tau_l K(s_i, s_{i'})$ with Gaussian kernel $K(s_i, s_{i'}) = \exp(-\frac{\|s_i - s_{i'}\|^2}{\gamma})$. $\tau_l$ and $\gamma$ are two hyperparameters controlling the scale of the factor value relative to the residual error and the strength of the spatial correlation, respectively. $\gamma$ is set numerically, in relation to the range of values in the data. $\tau_l = \tau$ is assumed constant across factors.

Maximum likelihood-based inference predicts the factor loading matrix $W$ and the factor matrix $Z$ together with the hyper-parameters $(\tau, \sigma_0^2)$.

The spatial domains are obtained by computing a clustering algorithm on $Z$, in this case Walktrap algorithm and Louvain algorithm for clustering analysis.

**Randomized Spatial PCA (RASP): a computationally efficient method for dimensionality reduction of high-resolution spatial transcriptomics data [32]** This approach first computes a random-
ized PCA on the expression matrix X with rank d and computes W, the matrix containing the top-d principal components of the covariance. The data is projected on the principal components and stored in a matrix Z = XW. In a second step, the spatial information is taken into account with the connectivity matrix C of the k-nearest neighbors graph. From this connectivity matrix, a distance matrix D is computed with diagonal elements set to a value α. The matrix is inverted with a power β, and stored in a matrix Dinv = (D⁻¹)β which is subsequently normalized column-wise. The final representation is obtained by computing the product Zs = DinvZ = (D⁻¹)βXW, which gives a spatially smoothed version of PCA-derived projection.

## Graph-based spatial PCA

**GraphPCA: a fast and interpretable dimension reduction algorithm for spatial transcriptomics data** [16] The expression matrix X is assumed to have been normalized by analytic Pearson residuals proposed by [45] and further scaled for each gene to have zero mean and unit standard deviation. Then, the top 3000 spatial variable genes of the SPARK package are selected.

GraphPCA takes physical distance into account by including a graph regularization to the variational PCA (3) leading to

$$ \min_{Z,W} \|X - ZW^T\|_2^2 + \frac{1}{2}\lambda \sum_{i,j \in C} C_{i,j} \|Z_i - Z_j\|_2^2, \text{ with } W^T W = I \quad (6) $$

where λ is a hyperparameter controlling the reconstruction error versus the smoothness of the projections on the underlying spatial neighborhood graph. This optimization problem can be rewritten as

$$ \min_{Z,W} \|X - ZW^T\|_2^2 + \lambda \mathrm{Tr}(Z^T \mathrm{Lap} Z), \text{ with } W^T W = I \quad (7) $$

where Lap = D - C is the Laplacian matrix of the spatial neighborhood graph, with D the diagonal matrix containing the degree of the nodes in the graph C.

The solution of equation (7) obtained with Lagrange multiplier is Z* = (I + λLap)^-1 XW and from this it follows that W* is the collection of top-d eigenvectors of the graph-revised covariance matrix X^T(I + λLap)^-1X = X^T(I + λ(D - C))^{-1}X.

**PCA-based spatial domain identification with state-of-the-art performance** [33] The preprocessing step for the expression data consists in normalization and log-transform per cell. In the same way as in jsPCA and graphPCA, for nichePCA [33], the spatial information is encoded in a k-nearest neighbor algorithm. The joint spatial and expression representation is obtained by computing the mean gene expression of a cell and its neighbors. Classical PCA is computed on this joint representation to obtained a spatially aware low-dimensional representation.

## Comparison with jsPCA

We have reviewed here recent approaches aimed at adapting classical PCA to spatially aware dimensionality reduction for spatial transcriptomics data. While they revolve around similar ideas, the mathematical formulations differ, leading to different numerical solutions.
jsPCA is based, for a single ST sample, on the diagonalization of the following spatial covariance which is obtained as the product of the spatial autocorrelation and the gene expression covariance:

$$A = \frac{1}{2n} X^T (L + L^T) X = \frac{1}{2n} X^T (D^{-1} C + (D^{-1} C)^T) X \quad (8)$$

The corresponding embedding is obtained as $Z = XW$, where $W$ correspond to the top-$d$ eigenvectors of the spatial covariance $A$.

From this analytical solution, we can directly see that this formulation is quite different from the one called spatialPCA in [15], equations (4) and (5). Compared to jsPCA, spatialPCA is distance-based and relies on a Gaussian kernel for the spatial covariance. Moreover, the factors of spatialPCA [15] can be found with maximum-likelihood based inference. Their complete derivation is quite heavy and can be found in [15] supplementary material. It is noticeable that the influence of the spatial information is controlled by the hyperparameter $\tau$.

Similarly for the method RASP [32], the optimal representation is obtained by $Z_s = (D^{-1})^\beta XW$, where $D$ is a $k$-N spatial distance matrix, and $W$ corresponds to the eigenvectors of the PCA. Analytically, the solution is quite different from jsPCA and rely on a parameter $\beta$ to control the influence of space, as well as several heuristic steps.

GraphPCA is closer to jsPCA as it is a graph-based approach, relying on a Laplacian regularization of the variational formulation of PCA [46]. The corresponding graph-revised covariance matrix reads $X^T(I + \lambda(D - C))^{-1}X$, to be compared with the spacial covariance of jsPCA in equation (8). With this rewriting of the two covariance matrices, it is clear that the two approaches have rather different grounds.

Finally, NichePCA computes the mean of the expression of spatial location and its direct neighbor as the new feature on which to apply PCA. Although it exploits the direct neighbors information, this approach does not capture the full extent of cross terms that jsPCA captures by computing the spatial autocorrelation on the $k$-NN graph.

## References

[1] Rao A, Barkley D, França GS, Yanai I. Exploring tissue architecture using spatial transcriptomics. Nature. 2021;596(7871):211-20.

[2] Liu Z, Sun D, Wang C. Evaluation of cell-cell interaction methods by integrating single-cell RNA sequencing data with spatial information. Genome Biology. 2022;23(1):218.

[3] Ruan Z, Zhou W, Liu H, Wei J, Pan Y, Yan C, et al. Precise detection of cell-type-specific domains in spatial transcriptomics. Cell Reports Methods. 2024;4(8).

[4] Zhou R, Yang G, Zhang Y, Wang Y. Spatial transcriptomics in development and disease. Molecular Biomedicine. 2023;4(1):32.

[5] Oliveira MFd, Romero JP, Chung M, Williams SR, Gottscho AD, Gupta A, et al. High-definition spatial transcriptomic profiling of immune cell populations in colorectal cancer. Nature Genetics. 2025:1-12.

[6] Sampath Kumar A, Tian L, Bolondi A, Hernández AA, Stickels R, Kretzmer H, et al. Spatiotemporal transcriptomic maps of whole mouse embryos at the onset of organogenesis. Nature Genetics. 2023;55(7):1176-85.

[7] Jin Y, Zuo Y, Li G, Liu W, Pan Y, Fan T, et al. Advances in spatial transcriptomics and its applications in cancer research. Molecular Cancer. 2024;23(1):129.
[8] Moses L, Pachter L. Museum of spatial transcriptomics. Nature methods. 2022;19(5):534-46.

[9] Yuan Z, Zhao F, Lin S, Zhao Y, Yao J, Cui Y, et al. Benchmarking spatial clustering methods with spatially resolved transcriptomics data. Nature Methods. 2024;21(4):712-22.

[10] Dong K, Zhang S. Deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder. Nature communications. 2022;13(1):1739.

[11] Ren H, Walker B, Cang Z, Nie Q. Identifying multicellular spatiotemporal organization of cells with SpaceFlow. Nat. Commun. 13, 4076; 2022.

[12] Li J, Chen S, Pan X, Yuan Y, Shen Hb. CCST: Cell clustering for spatial transcriptomics data with graph neural network. 2021.

[13] Li Z, Zhou X. BASS: multi-scale and multi-sample analysis enables accurate cell type clustering and spatial domain detection in spatial transcriptomic studies. Genome biology. 2022;23(1):168.

[14] Zhao E, Stone MR, Ren X, Guenthoer J, Smythe KS, Pulliam T, et al. Spatial transcriptomics at subspot resolution with BayesSpace. Nature biotechnology. 2021;39(11):1375-84.

[15] Shang L, Zhou X. Spatially aware dimension reduction for spatial transcriptomics. Nature communications. 2022;13(1):7203.

[16] Yang J, Wang L, Liu L, Zheng X. GraphPCA: a fast and interpretable dimension reduction algorithm for spatial transcriptomics data. Genome Biology. 2024;25(1):287.

[17] Chen C, Kim HJ, Yang P. Evaluating spatially variable gene detection methods for spatial transcriptomics data. Genome Biology. 2024;25(1):18.

[18] Svensson V, Teichmann SA, Stegle O. SpatialDE: identification of spatially variable genes. Nature methods. 2018;15(5):343-6.

[19] Sun S, Zhu J, Zhou X. Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies. Nature methods. 2020;17(2):193-200.

[20] Zhu J, Sun S, Zhou X. SPARK-X: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies. Genome biology. 2021;22(1):184.

[21] Hao M, Hua K, Zhang X. SOMDE: a scalable method for identifying spatially variable genes with self-organizing map. Bioinformatics. 2021;37(23):4392-8.

[22] Dries R, Zhu Q, Dong R, Eng CHL, Li H, Liu K, et al. Giotto: a toolbox for integrative analysis and visualization of spatial expression data. Genome biology. 2021;22(1):78.

[23] Weber LM, Saha A, Datta A, Hansen KD, Hicks SC. nnSVG for the scalable identification of spatially variable genes using nearest-neighbor Gaussian processes. Nature communications. 2023;14(1):4059.

[24] Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573-87.

[25] Jombart T, Devillard S, Dufour AB, Pontier D. Revealing cryptic spatial patterns in genetic variability by a new multivariate method. Heredity. 2008;101(1):92-103.

[26] Guijarro C, Song S, Aigouy B, Clément R, Villoutreix P, Kelly RG. Single-cell morphometrics reveals T-box gene-dependent patterns of epithelial tension in the Second Heart field. Nature Communications. 2024;15(1):9512.

[27] Theis FJ, Cason TP, Absil PA. Soft dimension reduction for ICA by joint diagonalization on the Stiefel manifold. In: International Conference on Independent Component Analysis and Signal Separation. Springer; 2009. p. 354-61.

[28] Townsend J, Koep N, Weichwald S. Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation. Journal of Machine Learning Research. 2016;17(137):1-5.

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

[30] Maulana I, Roestam R. Optimizing KNN Algorithm Using Elbow Method for Predicting Voter Participation Using Fixed Voter List Data (DPT). Jurnal Sosial Teknologi. 2024;4(7):441-51.
[31] Montano V, Jombart T. An Eigenvalue test for spatial principal component analysis. BMC bioinformatics. 2017;18:1-7.

[32] Gingerich IK, Goods BA, Frost HR. Randomized Spatial PCA (RASP): a computationally efficient method for dimensionality reduction of high-resolution spatial transcriptomics data. Research Square. 2025:rs-3.

[33] Schaub DP, Yousefi B, Kaiser N, Khatri R, Puelles VG, Krebs CF, et al. PCA-based spatial domain identification with state-of-the-art performance. Bioinformatics. 2025;41(1):btaf005.

[34] Hastie T, Tibshirani R, Friedman J, et al. The elements of statistical learning. Springer series in statistics New-York; 2009.

[35] Hu J, Li X, Coleman K, Schroeder A, Ma N, Irwin DJ, et al. SpaGCN: Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network. Nature methods. 2021;18(11):1342-51.

[36] Yan G, Hua SH, Li JJ. Categorization of 34 computational methods to detect spatially variable genes from spatially resolved transcriptomics data. Nature Communications. 2025;16(1):1141.

[37] Maynard KR, Collado-Torres L, Weber LM, Uytingco C, Barry BK, Williams SR, et al. Transcriptome-scale spatial gene expression in the human dorsolateral prefrontal cortex. Nature neuroscience. 2021;24(3):425-36.

[38] Traag VA, Waltman L, Van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. Scientific reports. 2019;9(1):1-12.

[39] Chen A, Liao S, Cheng M, Ma K, Wu L, Lai Y, et al. Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays. Cell. 2022;185(10):1777-92.

[40] Haghverdi L, Lun AT, Morgan MD, Marioni JC. Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nature biotechnology. 2018;36(5):421-7.

[41] Vandereyken K, Sifrim A, Thienpont B, Voet T. Methods and applications for single-cell and spatial multi-omics. Nature Reviews Genetics. 2023;24(8):494-515.

[42] Cardoso JF, Souloumiac A. Jacobi angles for simultaneous diagonalization. SIAM journal on matrix analysis and applications. 1996;17(1):161-4.

[43] Lehoucq RB, Sorensen DC, Yang C. ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. SIAM; 1998.

[44] Congedo M, Phlypo R, Pham DT. Approximate joint singular value decomposition of an asymmetric rectangular matrix set. IEEE Transactions on Signal Processing. 2010;59(1):415-24.

[45] Lause J, Berens P, Kobak D. Analytic Pearson residuals for normalization of single-cell RNA-seq UMI data. Genome biology. 2021;22(1):258.

[46] Belkin M, Niyogi P. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation. 2003;15(6):1373-96.