JGR-NMF: joint graph-regularized non-negative matrix factorization for spatial domain identification

Overview of JGR-NMF

The JGR-NMF method integrates gene expression profiles with spatial coordinates to analyze spatial transcriptomics data. As illustrated in Fig. 1, the analysis begins with data preprocessing, where genes with low expression are filtered, read counts are normalized, and highly variable genes are selected. Subsequently, an adaptive gene expression similarity graph and a fixed spatial proximity graph are constructed to model relationships between spots. These graphs are incorporated into a joint non-negative matrix factorization (NMF) framework, enabling the expression matrix to be factorized while preserving both transcriptional similarity and spatial neighborhood structures. Finally, the low-dimensional embeddings are clustered to identify spatially coherent tissue domains. By jointly leveraging molecular and spatial constraints, JGR-NMF demonstrates improved performance in tissue domain discovery.

Data description

Four spatial transcriptomics datasets were used in this study, including two cancer datasets (Breast Cancer-1 dataset and Breast Cancer-2 dataset), one Mouse Kidney dataset and one Mouse Embryo dataset. The Breast Cancer-1 dataset was obtained from the 10x Visium platform, with ground truth labels based on the 20 regions defined by Xu et al. (2024). The Breast Cancer-2 dataset was also derived from the 10x Visium platform, primarily consisting of ductal carcinoma in situ and invasive carcinoma, with its ground truth labels divided into eight categories using the platform-provided GraphClust algorithm. The Mouse Kidney dataset with six ground truth labels was also downloaded from the 10x Visium platform. The Mouse Embryo dataset was provided by Lohoff et al. (2022). Unlike the previous three datasets, the Mouse Embryo data was generated by the seqFISH platform, offering single-cell resolution and containing 23 cell types. The basic information of the four datasets is presented in Table 1.

Data preprocessing

To ensure data quality and provide a reliable foundation for downstream analyses, spatial transcriptomics data are preprocessed through four main steps: gene filtering, data normalization, highly variable gene (HVG) selection, and data standardization.

Let X = [X₁, X₂, …, X_i, …, X_s] ∈ ℝ^m×s denote the gene expression count matrix, where m is the number of genes and s is the number of spatial spots. Each column vector X_i = [X_i1, X_i2, …, X_im]^T represents the gene expression profile of the i-th spot. Gene filtering is first applied to remove low-quality or uninformative genes. A gene is discarded if it is not expressed in more than 80% of the spots. The zero-expression ratio for gene j is computed as $ze{x}_j=\frac{1}{s}{\sum }_{i=1}^s\mathbb{I}\left(X_{ij}=0\right)$, where 𝕀(⋅) is the indicator function that returns 1 if X_ij = 0, and 0 otherwise. Genes with zex_j ≥ τ (threshold τ = 0.8) are filtered out, resulting in a reduced matrix X′ ∈ ℝ^m′×s.

To eliminate technical variability between spots, the filtered expression matrix is normalized so that the gene expression values in each spot sum to one: $X_{ij}^{\prime }\leftarrow X_{ij}^{\prime }/{\sum }_{j=1}^{m^{\prime }}{X}_{ij}^{\prime }$. The resulting matrix X′ is used for HVG selection. Highly variable genes are selected to reduce dimensionality while preserving biological signal. For each gene j, the expression variance is calculated as ${\sigma }_j^2=\frac{1}{s}{\sum }_{i=1}^s{\left(X_{ij}^{\prime }-{\mu }_j\right)}^2$, where ${\mu }_j=\frac{1}{s}{\sum }_{i=1}^s{X}_{ij}^{\prime }$ is the mean expression of gene j. The top $\tilde{m}$ genes with the highest variances (typically $\tilde{m}=2000$) are retained, resulting in the HVG matrix $\tilde{X}\in {\mathbb{R}}^{\tilde{m}\times s}$.

Finally, the HVG matrix $\tilde{X}$ is standardized by scaling the expression profile of each gene across spots to the range [0, 1], using the transformation ${\tilde{X}}_{ij}\leftarrow {\tilde{X}}_{ij}/max\left({\tilde{X}}_j\right)$, where $max\left({\tilde{X}}_j\right)$ is the maximum expression value of gene j across all spots. The resulting standardized matrix $\tilde{X}$ is ready for downstream analyses.

Adaptive neighborhood graph based on spot relationships

The gene expression-based neighborhood graph in spatial transcriptomics data helps identify spots with similar gene expression patterns. Traditional neighborhood graphs typically require each data point to connect to a fixed number of nearest neighbors, regardless of whether those neighbors are located at the center or boundary of a cluster. However, for data points near the cluster boundary, their nearest neighbors may belong to unrelated clusters. Inspired by Cai, Huang & Yin (2022) and Li, Shyr & Liu (2024), an adaptive neighborhood graph is constructed that automatically adjusts the number of neighbors for each spot.

To characterize the positional relationships between spots, the gene expression distance between spots i and i′ is defined as $d_{i{i}^{\prime }}^1=\parallel {X^{\prime }}_i-{X^{\prime }}_{i^{\prime }}\parallel =\sqrt{{\sum }_{j=1}^m{\left({X^{\prime }}_{ij}-{X^{\prime }}_{i^{\prime }j}\right)}^2}$. Based on this expression distance, the neighbor set of spot i can be derived.

Let E = (E_ii′)_s×s = [E₁, E₂, …, E_i, …, E_s] denote the probability matrix representing neighbor connections between spots, where E_i = [E_i1, …, E_ii′, …, E_is]^T is the connection probability vector from spot i to all other spots. The value E_ii′ represents the probability that spot i connects to spot i′, and is inversely proportional to $d_{i{i}^{\prime }}^1$—the smaller $d_{i{i}^{\prime }}^1$, the larger E_ii′.

A maximum number of nearest neighbors $k_{max}^1$ is set for all spots to compute the distance distribution within the maximum neighborhood and adjust the actual number of neighbors according to sudden changes in distance. If the distance to the $\left(k_i^1+1\right)$-th neighbor is significantly larger than that to the $k_i^1$-th neighbor, then the true number of neighbors for spot i is $k_i^1$ ($k_i^1<k_{max}^1$). Figure 2 shows an example of adaptive k-nearest neighbors. In practice, the connection probability is employed instead of distance to identify the true number of nearest neighbors.

Following Cai, Huang & Yin (2022), the probabilities of spots connecting to their neighbors are computed by minimizing an objective function defined as follows. (1)${\displaystyle \underset{E_{i{i}^{\prime }}}{min}\sum _{i=1}^s\sum _{i^{\prime }=1}^s\sum _{n=0}^{\infty }\left(d_{i{i}^{\prime }}^1{E}_{i{i}^{\prime }}^n\right)\text{s.t.}0\le E_{i{i}^{\prime }}\le 1,\sum _{i^{\prime }=1}^s{E}_{i{i}^{\prime }}=1,i=1,2,\dots ,s.}$ Since ${\sum }_{n=0}^{\infty }{E}_{i{i}^{\prime }}^n={lim}_{n\to \infty }\frac{1-E_{i{i}^{\prime }}^n}{1-E_{i{i}^{\prime }}}\approx \frac{1}{1-E_{i{i}^{\prime }}}$ (Wu et al., 2019), Eq. (1) simplifies to (2)${\displaystyle \begin{array}{cc}\underset{E_{i{i}^{\prime }}}{min}\hfill & \sum _{i=1}^s\sum _{i^{\prime }=1}^s\frac{d_{i{i}^{\prime }}^1}{1-E_{i{i}^{\prime }}}\hfill \\ \text{s.t.}\hfill & 0\le E_{i{i}^{\prime }}\le 1,\hfill \\ \hfill & \sum _{i^{\prime }=1}^s{E}_{i{i}^{\prime }}=1,i=1,2,\dots ,s.\hfill \end{array}}$ The Lagrangian function of the minimization problem (2) is given by ${\displaystyle L\left(E,\zeta ,\lambda \right)=\sum _{i=1}^s\sum _{i^{\prime }=1}^s\frac{d_{i{i}^{\prime }}^1}{1-E_{i{i}^{\prime }}}-\sum _{i=1}^s{\zeta }_i\left(\sum _{i^{\prime }=1}^s{E}_{i{i}^{\prime }}-1\right)-\sum _{i=1}^s\sum _{i^{\prime }=1}^s{\lambda }_{i{i}^{\prime }}{E}_{i{i}^{\prime }},}$ where ζ and λ are Lagrange multipliers. Taking the derivative of L(E, ζ, λ) with respect to E_ii′ and setting it to zero, we get ${\displaystyle \frac{\partial L}{\partial E_{i{i}^{\prime }}}=\frac{d_{i{i}^{\prime }}^1}{{\left(1-E_{i{i}^{\prime }}\right)}^2}-{\zeta }_i-{\lambda }_{i{i}^{\prime }}=0.}$ Using the Karush-Kuhn–Tucker (KKT) conditions (Chen, Zeng & Pan, 2022), the optimal solution to the problem is given by $E_{i{i}^{\prime }}={\left(1-\frac{\sqrt{d_{i{i}^{\prime }}^1}}{\sqrt{{\zeta }_i}}\right)}_{+}$, where (⋅)₊ denotes the positive part, i.e., values where E_ii′ > 0. According to the imposed constraints, it follows that ${\displaystyle \sum _{i^{\prime }=1}^s{\left(1-\frac{\sqrt{d_{i{i}^{\prime }}^1}}{\sqrt{{\zeta }_i}}\right)}_{+}=1,\text{for}i=1,2,\dots ,s.}$

Assume that spot i has $k_i^1$ nearest neighbors such that $\left(1-\frac{\sqrt{d_{i{i}^{\prime }}^1}}{\sqrt{{\zeta }_i}}\right)>0$, and each spot can have at most $k_{\text{max}}^1$ neighbors. Let the distances from spot i to all other spots be sorted in ascending order as $\left[{\widehat{d}}_{i1}^1,\dots ,{\widehat{d}}_{il}^1,\dots ,{\widehat{d}}_{is}^1\right]$. Then, the constraint can be rewritten as ${\sum }_{l=1}^{k_i^1}\left(1-\frac{\sqrt{{\widehat{d}}_{il}^1}}{\sqrt{{\zeta }_i}}\right)=1.$ Define ${\sum }_{l=k_i^1+1}^{k_{\text{max}}^1}\left(1-\frac{\sqrt{{\widehat{d}}_{il}^1}}{\sqrt{{\zeta }_i}}\right)=\sigma ,$ where −1 < σ ≤ 0. Then we have ${\sum }_{l=1}^{k_{\text{max}}^1}\left(1-\frac{\sqrt{{\widehat{d}}_{il}^1}}{\sqrt{{\zeta }_i}}\right)=1+\sigma .$

Solving for $\sqrt{{\zeta }_i}$ yields $\sqrt{{\zeta }_i}=\frac{{\sum }_{l=1}^{k_{\text{max}}^1}\sqrt{{\widehat{d}}_{il}^1}}{k_{\text{max}}^1-1-\sigma }.$ Substituting this expression into the optimal solution, we obtain ${\displaystyle E_{i{i}^{\prime }}=\left\{\begin{array}{cc}\frac{\sum _{l=1}^{k_{\text{max}}^1}\sqrt{{\widehat{d}}_{il}^1}-\left(k_{\text{max}}^1-1-\sigma \right)\sqrt{d_{i{i}^{\prime }}^1}}{\sum _{l=1}^{k_{\text{max}}^1}{\widehat{d}}_{il}^1},\hfill & \text{if}{d}_{i{i}^{\prime }}^1\le {\widehat{d}}_{i{k}_{\text{max}}^1}^1,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\right.}$

Here, ${\widehat{d}}_{i{k}_{max}^1}^1$ represents the distance between spot i and its $k_{max}^1$-th nearest neighbor. The condition $d_{i{i}^{\prime }}^1\le {\widehat{d}}_{i{k}_{max}^1}^1$ ensures that only the top $k_{max}^1$ neighbors receive nonzero connection probabilities. Since E_ii′ > 0, the number of selected neighbors $k_i^1$ for spot i can be determined from Li, Shyr & Liu (2024) as (3)${\displaystyle k_i^1=\left\{\begin{array}{cc}{k}_{\text{max}}^1,\hfill & {\widehat{d}}_{k_{\text{max}}^1}^1<{\left(\frac{\sum _{l=1}^{k_{\text{max}}}\sqrt{{\widehat{d}}_{il}^1}}{k_{\text{max}}^1-1-\sigma }\right)}^2,\hfill \\ l^{\prime },\hfill & {\widehat{d}}_{i{l}^{\prime }}^1<{\left(\frac{\sum _{l=1}^{k_{\text{max}}}\sqrt{{\widehat{d}}_{il}^1}}{k_{\text{max}}^1-1-\sigma }\right)}^2\text{and}{\widehat{d}}_{i\left(l^{\prime }+1\right)}^1<{\left(\frac{\sum _{l=1}^{k_{\text{max}}}\sqrt{{\widehat{d}}_{il}^1}}{k_{\text{max}}^1-1-\sigma }\right)}^2,\hfill \end{array}\right.}$ where $l^{\prime }<k_{\text{max}}^1$. Based on Eq. (3), we obtain the neighbor set $N_i^1$ for spot i. The adjacency matrix $A^1={\left(A_{i{i}^{\prime }}^1\right)}_{s\times s}$ of the adaptive neighborhood graph based on gene expression is defined as follows (4)${\displaystyle A_{i{i}^{\prime }}^1=\left\{\begin{array}{cc}\frac{\sum _{l=1}^{k_{max}}\sqrt{{\widehat{d}}_{il}^1}-\left(k_{max}^1-1-\sigma \right)\sqrt{d_{i{i}^{\prime }}^1}}{\sum _{l=1}^{k_{max}^1}{\widehat{d}}_{il}^1},\hfill & \text{if}{i}^{\prime }\in N_i^1\text{and}i\in N_{i^{\prime }}^1,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\right.}$

The value of σ is selected via grid search. Specifically, σ is set to different values (0.1, 0.2, …, 0.9), and Louvain clustering is performed based on the resulting A¹. The optimal σ is chosen as the one just before a sharp increase in the number of identified communities.

Joint graph-regularized nonnegative matrix factorization model

Let Y = [Y₁, Y₂, …, Y_i, …, Y_s] ∈ ℝ^2×s denote the spatial coordinate matrix corresponding to the spatial transcriptomics data X, where Y_i = [Y_i1, Y_i2]^T represents the 2D coordinate vector of spot i. The spatial distance between spots i and i′ is defined as $d_{i{i}^{\prime }}^2$, given by $d_{i{i}^{\prime }}^2=\parallel Y_i-Y_{i^{\prime }}\parallel =\sqrt{{\left(Y_{i1}-Y_{i^{\prime }1}\right)}^2+{\left(Y_{i2}-Y_{i^{\prime }2}\right)}^2}.$ The set of k² nearest neighbors of spot i based on spatial location, denoted $N_i^2$, can be obtained. Define $A^2={\left(A_{i{i}^{\prime }}^2\right)}_{s\times s}$ as the adjacency matrix for spatial positions, where (5)${\displaystyle A_{i{i}^{\prime }}^2=\left\{\begin{array}{cc}1,\hfill & \text{if}{i}^{\prime }\in N_i^2\text{and}i\in N_{i^{\prime }}^2,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\right.}$

Let $L^1={\left(L_{i{i}^{\prime }}^1\right)}_{s\times s}$ denote the graph Laplacian matrix that reflects the topological relationships between gene expression neighbors. Specifically, (6)${\displaystyle L^1=D^1-A^1,}$ (7)${\displaystyle D_{i{i}^{\prime }}^1=\left\{\begin{array}{cc}\sum _{i^{\prime }}{A}_{i{i}^{\prime }}^1,\hfill & i=i^{\prime },\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.}$

Similarly, let $L^2={\left(L_{i{i}^{\prime }}^2\right)}_{s\times s}$ be the Laplacian matrix describing the topological relationships among spatial neighbors (8)${\displaystyle L^2=D^2-A^2,}$ (9)${\displaystyle D_{i{i}^{\prime }}^2=\left\{\begin{array}{cc}\sum _{i^{\prime }}{A}_{i{i}^{\prime }}^2,\hfill & i=i^{\prime },\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.}$

To obtain a low-dimensional representation of the spatial transcriptomics data, we first process the gene expression matrix X′ by selecting $\tilde{m}$ highly variable genes (default $\tilde{m}=2000$), and then normalize the matrix of selected genes to obtain $\tilde{X}\in {\mathbb{R}}^{\tilde{m}\times s}$. Based on this, we propose the Joint Graph-Regularized Nonnegative Matrix Factorization (JGR-NMF) model that integrates gene expression and spatial information (10)${\displaystyle \underset{W\ge 0,H\ge 0}{min}{\parallel \tilde{X}-WH\parallel }_F^2+{\alpha }_1\text{Tr}\left(H{L}^1{H}^T\right)+{\alpha }_2\text{Tr}\left(H{L}^2{H}^T\right),}$ where ${\parallel \cdot \parallel }_F^2$ denotes the Frobenius norm, and α₁ and α₂ are positive regularization parameters. $W\in {\mathbb{R}}^{\tilde{m}\times f}$ and H ∈ ℝ^f×s are the basis matrix and coefficient matrix, respectively.

This model incorporates two adjacency matrices—A¹ (adaptive adjacency matrix based on gene expression) and A² (adjacency matrix based on spatial coordinates)—to simultaneously capture expression similarity and spatial proximity, thereby providing a more comprehensive representation of latent structures and spatial heterogeneity in spatial transcriptomics data. Notably, the matrix H represents the low-dimensional embedding of the spatial transcriptomics data. Therefore, clustering can be performed by applying the K-means algorithm to H.

Model optimization

Following previous studies (Qiu et al., 2024; Wang, Bensmail & Gao, 2013), we adopt an alternating update algorithm to solve the optimization problem Eq. (10). First, we construct the following Lagrangian function (11)${\displaystyle \begin{array}{cc}\mathcal{L}\hfill & ={\parallel \tilde{X}-WH\parallel }_F^2+{\alpha }_1\text{Tr}\left(H{L}^1{H}^T\right)+{\alpha }_2\text{Tr}\left(H{L}^2{H}^T\right)+\text{Tr}\left(\Psi W^T\right)+\text{Tr}\left(\Phi H^T\right)\hfill \\ \hfill & =\text{Tr}\left(\tilde{X}{\tilde{X}}^T\right)-2\text{Tr}\left(H^T{W}^T\tilde{X}\right)+\text{Tr}\left(WH{H}^T{W}^T\right)+{\alpha }_1\text{Tr}\left(H{L}^1{H}^T\right)\hfill \\ \hfill & +{\alpha }_2\text{Tr}\left(H{L}^2{H}^T\right)+\text{Tr}\left(\Psi W^T\right)+\text{Tr}\left(\Phi H^T\right),\hfill \end{array}}$ where Ψ and Φ are the Lagrange multipliers associated with the constraints W ≥ 0 and H ≥ 0, respectively. Taking the partial derivatives of the Lagrangian function $\mathcal{L}$ with respect to W and H, we obtain (12)${\displaystyle \frac{\partial \mathcal{L}}{\partial W}=-2\tilde{X}{H}^T+2WH{H}^T+\Psi ,}$ (13)${\displaystyle \frac{\partial \mathcal{L}}{\partial H}=-2W^T\tilde{X}+2W^{T}WH+2H\left({\alpha }_1{L}^1+{\alpha }_2{L}^2\right)+\Phi .}$

Setting Eqs. (12) and (13) to zero and multiplying both sides by W^T and H^T, respectively, yield (14)${\displaystyle -2\tilde{X}{H}^T{W}^T+2WH{H}^T{W}^T+\Psi W^T=0,}$ (15)${\displaystyle -2W^T\tilde{X}{H}^T+2W^{T}WH{H}^T+2H\left({\alpha }_1{L}^1+{\alpha }_2{L}^2\right)H^T+\Phi H^T=0.}$

Substituting Eqs. (6) and (8) into (15) gives (16)${\displaystyle -2\left(W^T\tilde{X}+{\alpha }_{1}H{A}^1+{\alpha }_{2}H{A}^2\right)H^T+2\left(W^{T}WH+{\alpha }_{1}H{D}^1+{\alpha }_{2}H{D}^2\right)H^T+\Phi H^T=0.}$

According to the KKT conditions (Chen, Zeng & Pan, 2022), we have Ψ_pqW_pq = 0 and Φ_pqH_pq = 0. Therefore, Eqs. (14) and (16) can be rewritten as (17)${\displaystyle {\left(-\tilde{X}{H}^T+WH{H}^T\right)}_{pq}{W}_{pq}=0,}$ (18)${\displaystyle {\left(-W^T\tilde{X}-{\alpha }_{1}H{A}^1-{\alpha }_{2}H{A}^2+W^{T}WH+{\alpha }_{1}H{D}^1+{\alpha }_{2}H{D}^2\right)}_{pq}{H}_{pq}=0.}$

From Eqs. (17) and (18), the update rules for the basis matrix W and the coefficient matrix H can be derived as (19)${\displaystyle W_{pq}\leftarrow W_{pq}\frac{{\left(\tilde{X}{H}^T\right)}_{pq}}{{\left(WH{H}^T\right)}_{pq}},}$ (20)${\displaystyle H_{pq}\leftarrow H_{pq}\frac{{\left(W^T\tilde{X}+{\alpha }_{1}H{A}^1+{\alpha }_{2}H{A}^2\right)}_{pq}}{{\left(W^{T}WH+{\alpha }_{1}H{D}^1+{\alpha }_{2}H{D}^2\right)}_{pq}}.}$

As in previous works (Cai et al., 2011; Long et al., 2013), the update rules for dual-graph regularized non-negative matrix factorization are guaranteed to converge. Algorithm 1  provides the pseudocode for solving the proposed JGR-NMF.

 
_______________________ 
Algorithm 1 Joint Graph-Regularized Non-negative Matrix Factorization_____ 
Input: Gene expression count matrix X ∈ Rm×s; 
         Spatial location matrix Y ∈ R2×s; 
         Model parameters α1, α2; 
         Integers k1max, k2, m′, f; 
         Convergence threshold ϵ; 
         Maximum number of iterations MNI; 
Output: W and H 
 1:  Perform preliminary gene filtering and data normalization on X to obtain 
     the filtered matrix X′ ∈ Rm′×s 
      ; 
  2:  Select the top ^ m highly variable genes from X′ to obtain the reduced matrix 
      ^ X ∈ R^m×s; 
  3:  Construct adjacency matrices A1 and A2 using equations (4) and (5); 
  4:  Construct matrices D1 and D2 using equations (7) and (9); 
  5:  Construct Laplacian matrices L1 and L2 using equations (6) and (8); 
  6:  Randomly initialize W = W0 ∈ Rm′×f 
       ; 
  7:  Randomly initialize H = H0 ∈ Rf×s; 
  8:  for t = 1,2,...,MNI do 
 9:     Update W using equation (19); 
10:     Update H using equation (20); 
11:     Compute the error errH = ∥Ht − Ht−1∥∞; 
12:     if errH < ϵ then 
13:         Terminate iteration; 
14:         return  W and H; 
15:     end if 
16:  end for 
17:  Return: W and H______________________________________________________________    

Computational complexity analysis

To evaluate the efficiency of JGR-NMF in processing large-scale data, we conducted a computational complexity analysis. The cost of the algorithm mainly comes from the alternating iterative updates of the base matrix W and the coefficient matrix H. Specifically, according to the update rule (19) (Equation (19)), the computational complexity of updating W is O(m × s × f). Similarly, according to the update rule (20) (Equation (20)), the complexity of updating H is also O(m × s × f). Theoretically, constructing the Laplacian matrices L¹ and L² has a complexity of O(s²). However, since the constructed adjacency matrices A¹ and A² are highly sparse (each spot has an average of only k neighbors and k ≪ s), they can be optimized to O(s × k) through sparse matrix operation techniques. Therefore, the overall complexity of a single iteration can be controlled at O(m × s × f). Assuming the number of iterations required for algorithm convergence is t, the total computational complexity of JGR-NMF is O(t × m × s × f).

Comparison methods and evaluation metrics

To evaluate the performance of JGR-NMF, we benchmarked it against five existing methods: one non-spatial clustering approach (Seurat; Stuart et al., 2019) and four spatial clustering methods (SpaGCN; Hu et al., 2021), BANKSY (Singhal et al., 2024), MNMST (Wang, Liu & Ma, 2024), and spaMMCL (Liang et al., 2024b). For all methods, the number of clusters was set identically using default parameter configurations specified in their respective implementations. Performance was evaluated using three quantitative metrics: Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), and Purity (PUR). The definitions of these metrics are as follows:

(1) ARI. ARI measures the consistency between the clustering results and the ground truth labels, adjusting for the chance grouping of elements. It ranges from −1 to 1, with values closer to 1 indicating higher clustering consistency. The ARI is defined as: ${\displaystyle ARI=\frac{RI-E\left(RI\right)}{max\left(RI\right)-E\left(RI\right)},}$ where RI is the Rand Index, which reflects the similarity between the clustering results and true labels, E(RI) is the expected value of RI under random labeling, and max(RI) is the maximum value of RI. The Rand Index is computed as $RI=\frac{a+d}{a+b+c+d},$ where a is the number of sample pairs that are in the same cluster and have the same label, b is the number of pairs that are in the same label but in different clusters, c is the number of pairs that are in the same cluster but have different labels, and d is the number of pairs that are in different clusters and have different labels.

(2) NMI. NMI measures the amount of shared information between the predicted clusters and the ground truth labels, normalized to eliminate the influence of differing numbers of clusters and labels. It ranges from 0 to 1, with higher values indicating better clustering performance. Let there be s data points, and let O = {O₁, O₂, …, O_Kt} denote the set of ground truth labels and $\widehat{O}=\left\{{\widehat{O}}_1,{\widehat{O}}_2,\dots ,{\widehat{O}}_{K_f}\right\}$ denote the set of predicted clusters, where K_t is the number of true classes and K_f is the number of predicted clusters. Then NMI is defined as ${\displaystyle NMI\left(\widehat{O},O\right)=\frac{2I\left(\widehat{O},O\right)}{H\left(\widehat{O}\right)+H\left(O\right)},}$ where $I\left(\widehat{O},O\right)={\sum }_{u=1}^{K_t}{\sum }_{v=1}^{K_f}P\left(u,v\right)log\frac{P\left(u,v\right)}{P\left(u\right)P\left(v\right)}$ is the mutual information between $\widehat{O}$ and O, with P(u, v) denoting the joint probability that a sample belongs to both class O_u and cluster ${\widehat{O}}_v$, and P(u), P(v) denoting the marginal probabilities of O_u and ${\widehat{O}}_v$, respectively. The entropies H(O) and $H\left(\widehat{O}\right)$ are given by $H\left(O\right)=-{\sum }_{u=1}^{K_t}P\left(O_u\right)logP\left(O_u\right)$ and $H\left(\widehat{O}\right)=-{\sum }_{v=1}^{K_f}P\left({\widehat{O}}_v\right)logP\left({\widehat{O}}_v\right)$, respectively.

(3) PUR. Purity measures the extent to which each cluster contains samples from a single ground truth class. It ranges from 0 to 1, with higher values indicating more homogeneous clusters. It is defined as: ${\displaystyle PUR=\frac{1}{s}\sum _{v=1}^{K_f}\underset{u}{max}|{\widehat{O}}_v\cap O_u|,}$ where s is the total number of samples, K_f is the number of predicted clusters, ${\widehat{O}}_v$ is the set of samples in the v-th predicted cluster, and O_u is the set of samples with the u-th true label. The term $|{\widehat{O}}_v\cap O_u|$ represents the number of samples in cluster ${\widehat{O}}_v$ that also belong to true class O_u, and ${max}_u|{\widehat{O}}_v\cap O_u|$ gives the number of samples in the majority class within cluster v.

Spatial domain identification on the Breast Cancer-1 dataset

To evaluate the spatial domain identification performance of JGR-NMF, we first conducted comparative experiments on the Breast Cancer-1 dataset against Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL. Figure 3A shows the spatial distribution of 20 manually annotated ground truth labels in the Breast Cancer-1 dataset. Figure 3B shows the ARI, NMI and PUR scores of different methods on the Breast Cancer-1 dataset in a more intuitive bar chart format. The x-axis represents different evaluation metrics, the y-axis represents the scores, and different colors correspond to different methods. Higher values indicate better clustering performance.

The experimental results demonstrate that JGR-NMF achieved the highest ARI, NMI, and PUR scores on the Breast Cancer-1 dataset, indicating superior performance in spatial domain identification. Specifically, JGR-NMF achieved an ARI of 0.6210, which is 0.1599, 0.0452, 0.2097, 0.0523, and 0.0511 higher than those of Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL, respectively. Its NMI value reached 0.6889, surpassing the compared methods by 0.0780, 0.0221, 0.0483, 0.0174, and 0.0235. Additionally, JGR-NMF achieved a PUR value of 0.7270, exceeding Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL by 0.0487, 0.0342, 0.0187, 0.0094, and 0.0382, respectively. On the Breast Cancer-1 dataset, the hyperparameters α₁, α₂, $k_{max}^1$ and σ in JGR-NMF are 0.8, 0.1, 100 and 0.3 respectively. To evaluate the statistical significance of the performance difference of different methods on the Breast Cancer-1 dataset, we performed paired t-tests to quantify the difference between JGR-NMF and each baseline method in ARI, NMI and PUR indicators (Fig. 3B). The significance levels are marked as ∗ (p < 0.05), ∗∗ (p < 0.01), and ∗∗∗ (p < 0.001). The statistical test results indicate that JGR-NMF exhibits significant advantages over most baseline methods in terms of ARI, NMI, and PUR indicators.

Figure 3C further visualizes the spatial domain identification results of all six methods. These visualizations reveal the capacity of each method to recover the true spatial structure. JGR-NMF demonstrates a high degree of consistency with the ground truth labels in multiple regions, especially in “IDC_2”, “IDC_4”, “DCIS/LCIS 1”, “DCIS/LCIS 4”, and “Tumor_edge_1”, accurately capturing the spatial characteristics of these regions. In contrast, although other methods show some effectiveness in certain areas, their overall performance is inferior. For example, BANKSY performs relatively well in the “Healthy_1” region but suffers from fragmented clustering in other areas and fails to capture the overall spatial distribution of the true labels. Similarly, while SpaGCN and MNMST perform moderately well in some regions, they still fall short in complex labeling areas compared to JGR-NMF.

To further verify the biological rationality of the spatial domain identified by JGR-NMF on the Breast Cancer-1 dataset, we conducted differential expression analysis on the predicted spatial domain. Specifically, we adopted the non-parametric Wilcoxon rank-sum test to screen the significantly highly expressed genes that were expressed in at least 25% of the spots in each spatial domain and whose expression levels was |log2FoldChange| ≥ 0.25 as candidate marker genes. By comparing the correspondence between predicted clusters and real spatial domain labels, it was found that predicted cluster 5 is highly consistent in spatial distribution with the real histological label “IDC4” domain. Based on this, we screened the top 20 marker genes with the highest |log2FoldChange| values from cluster 5, and performed Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The GO analysis results in Fig. 4A show that these genes are significantly enriched in multiple biological processes closely related to tumor progression, including “peptide hormone processing”, “cell–cell adhesion via plasma-membrane adhesion molecules”, and “heart valve development”. Cellular component terms were primarily associated with extracellular matrix structures and motor protein complexes, such as “collagen-containing extracellular matrix” and “dynein complex”. Molecular functions included “minus-end-directed microtubule motor activity” and “serine-type endopeptidase activity”. The KEGG pathway analysis in Fig. 4B further shows that these genes are significantly enriched in “Cell adhesion molecules”, “TGF-beta signaling pathway”, “Pathways of neurodegeneration-multiple diseases” and other pathways closely related to the occurrence and development of breast cancer. These results not only confirm the effectiveness of JGR-NMF in identifying biologically consistent spatial domains, but also reveal the potential molecular mechanisms of the IDC4 domain and its functional role in the tumor microenvironment.

Spatial domain identification on the Breast Cancer-2 dataset

To evaluate the performance of the JGR-NMF method in identifying spatial domains under different clustering numbers, we further conducted experiments on the Breast Cancer-2 dataset. Figures 5A and 5B show the manually annotated reference regions and the eight spatial domains divided by the GraphClust algorithm based on the 10x Visium platform, respectively. Since the manually annotated regions do not perfectly match the spot labels in the data, the eight regions identified by GraphClust are used as ground truth labels to evaluate the performance of different spatial domain identification methods.

Figure 5C visually displays the ARI, NMI, and PUR values of different methods on the Breast Cancer-2 dataset through a bar chart. The experimental results show that JGR-NMF performs best on the spatial domain identification task for the Breast Cancer-2 dataset, significantly outperforming other methods in terms of ARI, NMI, and PUR. Specifically, JGR-NMF achieves an ARI of 0.5354, which is higher than Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL by 0.2479, 0.1876, 0.1334, 0.3163, and 0.3646, respectively. Its NMI value of 0.5599 exceeds those of Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL by 0.0888, 0.0725, 0.0582, 0.2006, and 0.3087, respectively. Likewise, the PUR of JGR-NMF is 0.7431, surpassing Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL by 0.1827, 0.1136, 0.1013, 0.2645, and 0.2780, respectively. For the Breast Cancer-2 dataset, we set the JGR-NMF hyperparameters to α₁ = 0.2, α₂ = 1, $k_{max}^1=50$, and σ = 0.5. In addition, the t-test results in Fig. 5C indicate that JGR-NMF is significantly superior to other comparative methods in terms of ARI, NMI, and PUR indicators.

The spatial domain identification visualizations of the six methods on the Breast Cancer-2 dataset are shown in Fig. 5D. From these visualizations, it is evident that the domains identified by JGR-NMF are more consistent with the ground truth labels, especially in fibrous tissue and invasive cancer regions corresponding to ground truth labels “2”, “3”, “4”, and “6” (Fig. 5A), where JGR-NMF exhibits strong identification ability for cancer-related areas.

In contrast, although some methods have partial overlap with the ground truth, their overall performance is suboptimal. For example, Seurat tends to produce overly concentrated clusters, with most spots assigned to a single cluster, leading to low discrimination. BANKSY, on the other hand, produces overly fragmented clusters and fails to form meaningful spatial structures. SpaGCN, MNMST, and spaMMCL tend to aggregate spatially adjacent regions into large blocks, resulting in boundaries that deviate from the ground truth.

In conclusion, JGR-NMF not only achieves the best performance in terms of quantitative evaluation metrics such as ARI, NMI, and PUR but also demonstrates superior spatial domain identification in visual inspection, outperforming Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL.

Spatial domain identification on the mouse kidney dataset

To further evaluate the generalization performance of JGR-NMF method in identifying spatial domains on different datasets, we conducted comparative experiments with Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL methods on the Mouse Kidney dataset. Figure 6A shows the spatial distribution of six ground truth labels manually annotated in the Mouse Kidney dataset. Similar to the visualization results on the breast cancer dataset, we used the bar chart in Fig. 6B to compare the ARI, NMI and PUR scores of different methods on the Mouse Kidney dataset. The experimental results indicate that JGR-NMF performs the best in the spatial domain identification task on the Mouse Kidney dataset. Specifically, the ARI value of JGR-NMF is 0.6074, which is 0.0333, 0.0718, 0.0624, 0.1416, and 0.1934 higher than Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL, respectively. The NMI value of JGR-NMF is 0.6463, which is 0.0124, 0.0345, 0.0399, 0.1078, and 0.1308 higher than Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL, respectively. Similarly, the PUR of JGR-NMF is 0.7531, which is 0.0222, 0.0174, 0.0083, and 0.1078 higher than SpaGCN, BANKSY, MNMST, and spaMMCL, respectively. Although Seurat has a slight advantage in purity indicators on the Mouse Kidney dataset, the JGR-NMF method is the best based on the comprehensive ARI and NMI indicators. On the Mouse Kidney dataset, the hyperparameters for JGR-NMF were set as follows: α₁ = 0.001, α₂ = 0.999, $k_{max}^1=45$, and σ = 0.6. In addition, the t-test results in Fig. 6B indicate that JGR-NMF is significantly superior to other comparative methods in ARI, NMI, and PUR indicators on the Mouse Kidney dataset.

The spatial domain identification results of all six methods are presented in Fig. 6C. It can be observed that JGR-NMF shows a high degree of consistency with the ground truth labels of multiple domains. In contrast, although other methods demonstrate some effectiveness in certain domains, their overall performance is poorer.

Spatial domain identification on the mouse embryo dataset

To evaluate the clustering performance of JGR-NMF on datasets from different spatial transcriptomics platforms, we conducted comparative experiments on the seqFISH Mouse Embryo dataset against Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL. The primary difference between seqFISH and 10x Visium data lies in the much higher single-cell resolution of seqFISH, which enables the direct detection of gene expression at the single-cell level while preserving spatial positional information. The high resolution of seqFISH data provides a suitable condition to verify the fine-grained performance of JGR-NMF.

Figure 7A shows the spatial distribution of 23 ground truth labels in the Mouse Embryo dataset. These regions exhibit curved, band-like, and relatively scattered spatial patterns. To better compare the performance of different methods, the ARI, NMI, and PUR scores of all methods are presented in Fig. 7B, where the horizontal axis represents different indicators and the vertical axis represents scores. Different methods are color-coded for clarity. The experimental results show that JGR-NMF achieves the highest PUR score on the Mouse Embryo dataset, and it also obtains the highest ARI and NMI values among all methods except for Seurat. Specifically, JGR-NMF achieves an ARI of 0.4478, which is higher than SpaGCN, BANKSY, MNMST, and spaMMCL by 0.0300, 0.0714, 0.3108, and 0.0589, respectively. In terms of NMI, JGR-NMF scores 0.6319, exceeding SpaGCN, BANKSY, MNMST, and spaMMCL by 0.0459, 0.0406, 0.3308, and 0.0457, respectively. Moreover, JGR-NMF achieves a PUR score of 0.7270, which is 0.0107, 0.0293, 0.0566, 0.3312, and 0.0674 higher than Seurat, SpaGCN, BANKSY, MNMST, and spaMMCL, respectively. On the Mouse Embryo dataset, we adopted the following hyperparameter configuration for JGR-NMF: α₁ = 0.1, α₂ = 0.6, $k_{max}^1=52$, and σ = 0.3. Paired t-test was adopted to evaluate the statistical significance of performance differences between different methods on a mouse embryo dataset. The statistical test results in Fig. 7B indicate that JGR-NMF exhibits significant advantages over most baseline methods in terms of ARI, NMI, and PUR indicators.

Furthermore, the spatial domain identification results of the six methods on the Mouse Embryo dataset are visualized in Fig. 7C, where different colors correspond to different labels. From the visual results, we observe that the spatial domains identified by JGR-NMF are more consistent with the ground truth labels and display clear band-like structures.

Although Seurat slightly outperforms JGR-NMF in ARI and NMI (by 0.0097 and 0.0179, respectively), this can be attributed to the fact that Seurat is a non-spatial clustering method and is inherently better suited for single-cell resolution data. Given that seqFISH data are essentially single-cell data, Seurat performs slightly better in this context.

However, it is worth noting that Fig. 7C shows that despite Seurat’s slight lead, its results do not effectively capture finer banded regions, whereas JGR-NMF performs better in this regard. Specifically, JGR-NMF clearly identifies multiple band-like regions that closely match the ground truth spatial layout. In contrast, MNMST shows significant clustering in red blood cell regions but produces blocky clusters in other areas, which deviate substantially from the ground truth. Additionally, SpaGCN, BANKSY, and spaMMCL exhibit weaker abilities in identifying band-like regions, with more scattered clustering results and poorer correspondence with the ground truth.

In conclusion, although Seurat shows a slight advantage in single-cell clustering tasks, JGR-NMF demonstrates stronger capabilities in identifying more complex spatial distribution patterns, especially in accurately recognizing spatial domains with complex banded structures.

The personal computer configuration used in the experiment is as follows: (1) CPU: Intel Core i5-8265U; (2) Memory: 16GB RAM; (3) Software: R (version 4.3.2; R Core Team, 2023) and Python (version 3.12). Tables 2 and 3 list the total runtime and memory usage of all methods on four datasets. On all datasets, the total running time of JGR-NMF is significantly shorter than that of deep learning based methods (SpaGCN and spaMMCL). Although the running time of JGR-NMF is slightly longer than Seurat, BANKSY, and SpaGCN, considering its significant advantage in clustering accuracy, this additional time cost is reasonable and acceptable in practical applications. In addition, JGR-NMF maintained a low memory usage in all experiments, especially on large-scale Mouse Embryo datasets. This was mainly due to the model’s approximately linear computational complexity and the optimized operations on sparse matrices.

Ablation study

To comprehensively evaluate the impact of two regularization terms on the performance of JGR-NMF, we conducted detailed comparative experiments. These regularization terms were constructed based on gene expression and spatial location information, respectively. Experiments were performed on multiple datasets, including the Breast Cancer-1 dataset, Breast Cancer-2 dataset, and the Mouse Embryo dataset. Let and denote the two regularization terms constructed according to Equations (6) and (8), respectively. By combining these regularization terms, we assessed the influence of different topological structures on the clustering performance of JGR-NMF.

Tables 4, 5 and 6 present the ablation study results on the Breast Cancer-1, Breast Cancer-2, and Mouse Embryo datasets, respectively, including three key metrics: ARI, NMI, and PUR. These metrics effectively reflect the clustering accuracy and spatial domain identification performance of different model variants. The experimental results show that the model integrating both gene expression-based and spatial position-based neighbor topologies—constructed using adaptive neighborhood mechanisms—achieves the highest ARI, NMI, and PUR scores across all datasets, demonstrating the advantage of the proposed method in handling complex spatial transcriptomics data.

On the Breast Cancer-1 dataset, the model with both regularization terms outperforms the model using only by 0.0068, 0.0076, and 0.0058 in ARI, NMI, and PUR, respectively, indicating that the integration of spatial and gene expression information significantly improves clustering results. This trend is further confirmed on the Breast Cancer-2 dataset, where ARI, NMI, and PUR increase by 0.0161, 0.0154, and 0.0327, respectively, showing that in more complex datasets, the integrated model can better exploit spatial structural features.

Moreover, on the Mouse Embryo dataset, the combined model also demonstrates superior performance, with ARI, NMI, and PUR increasing by 0.0056, 0.0085, and 0.0012, respectively, compared to using only . These results further verify the ability of model to stably improve clustering accuracy under different experimental settings. Overall, although the model constructed using only spatial information performs slightly better than using only gene expression, their combination significantly enhances the overall performance. This finding suggests that the joint use of spatial and gene expression information can effectively improve the applicability and accuracy of clustering methods on complex datasets, especially in tasks that require leveraging spatial distribution patterns, where the regularization-integrated model exhibits strong advantages.