Spatial transcriptomics iterative hierarchical clustering (stIHC): A novel method for identifying spatial gene co‐expression modules
Catherine Higgins, Jingyi Jessica Li, Michelle Carey
TL;DR
This paper introduces stIHC, a new method for clustering genes in spatial transcriptomics data to better identify spatial gene expression patterns and their biological functions.
Contribution
The novel stIHC method improves detection of spatial gene co-expression modules compared to existing approaches.
Findings
stIHC outperforms existing methods like SPARK and SpatialDE in clustering spatially variable genes.
Gene ontology analysis shows that stIHC modules are biologically relevant and functionally consistent.
The method is robust across different spatial transcriptomics technologies and resolutions.
Abstract
Recent advancements in spatial transcriptomics (ST) technologies allow researchers to simultaneously measure RNA expression levels for hundreds to thousands of genes while preserving spatial information within tissues, providing critical insights into spatial gene expression patterns, tissue organization, and gene functionality. However, existing methods for clustering spatially variable genes (SVGs) into co‐expression modules often fail to detect rare or unique spatial expression patterns. To address this, we present spatial transcriptomics iterative hierarchical clustering (stIHC), a novel method for clustering SVGs into co‐expression modules, representing groups of genes with shared spatial expression patterns. Through three simulations and applications to ST datasets from technologies such as 10x Visium, 10x Xenium, and Spatial Transcriptomics, stIHC outperforms clustering…
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7
FIGURE 8
FIGURE 9
FIGURE 10| Metric | stIHC | SPARK | SPARK‐X | MERINGUE | SpatialDE |
|---|---|---|---|---|---|
| ARI |
| 0.94 | 0.94 | 0.94 | 0.11 |
| DBI |
| 0.54 | 0.54 | 0.54 | 4.54 |
| Metric | stIHC | SPARK | SPARK‐X | MERINGUE | SpatialDE |
|---|---|---|---|---|---|
| ARI |
| 0.94 | 0.94 | 0.94 | 0.05 |
| DBI |
| 0.54 | 0.54 | 0.54 | 4.45 |
- —Science Foundation Ireland10.13039/501100001602
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsSingle-cell and spatial transcriptomics
INTRODUCTION
1
Recent advancements in spatial transcriptomics (ST) technologies have enabled the measurement of gene expression levels while preserving spatial information within tissues. These technologies include sequencing‐based platforms such as Spatial Transcriptomics [1], 10x Visium [2], and Slide‐seq [3, 4], as well as imaging‐based platforms such as 10x Xenium [5], MERFISH [6], and seqFISH [7], which allow the quantification of hundreds to thousands of genes in relation to their spatial location.
Many methods have been developed to identify spatially variable genes (SVGs), which exhibit significant spatial expression variability across a tissue. For a comprehensive overview of these methods, see Ref. [8]. Despite substantial progress in identifying SVGs, there remains a lack of methods for clustering these genes into co‐expression modules—groups of genes with similar spatial expression patterns. Since biological pathways and processes are often driven by the coordinated activity of multiple genes [9, 10, 11], identifying co‐expression modules can enhance our understanding of the spatial organization and function of tissues. For example, this approach can help identify regions in tumors with unique gene expression patterns and functional roles [12].
Numerous studies have focused on clustering cells or spatial spots in ST data; for an overview of these methods, see Ref. [13]. However, these methods differ fundamentally from our proposed method. Spot‐based methods cluster spatial locations, typically corresponding to cell types, states, or functional regions within the tissue. In contrast, our focus is on gene co‐expression modules, clustering at the gene level. Here, clusters represent groups of genes exhibiting similar spatial expression patterns across all spatial spots. Figure 1 illustrates this distinction. The left panel depicts the tissue slice, the middle panel overlays a grid representing the spatial spots where gene abundance is measured using ST technologies, and the right panel shows the resulting gene expression measurements for nine genes across the tissue slice. As shown, many genes exhibit similar spatial expression patterns; in this example, the nine genes display three distinct patterns. The aim of our method is to cluster genes with similar spatial behavior across the entire tissue into meaningful groups. This shift from spot‐based analysis to gene‐based analysis emphasizes uncovering relationships among genes rather than identifying cell types or spatial domains.
Workflow of spatial transcriptomics analysis. (Left) Tissue section, adapted from the Allen Mouse Brain Atlas (mouse.brain‐map.org) [14]. (Middle) Spatial transcriptomics grid overlay, representing spatially resolved gene expression measurements. (Right) Identified gene co‐expression modules, illustrating groups of genes with distinct spatial expression patterns. Each panel displays the spatial expression of a single gene.
Some existing methods for detecting SVGs in ST data incorporate gene clustering as a secondary step, but these methods primarily aim to identify SVGs rather than explicitly detect spatial gene co‐expression modules. Typically, these methods identify SVGs using a statistic calculated for each gene and subsequently group the detected SVGs into clusters. For example, SpatialDE [15] and SPARK [16] use Gaussian process (GP) regression to decompose the expression levels of each gene across a tissue slice into a spatial component and an independent noise term, identifying a gene as an SVG if the spatial component is statistically significant. The effectiveness of GP regression depends heavily on how accurately the selected spatial covariance matrix (derived from a chosen kernel) models spatial patterns [8]. After detecting SVGs, SpatialDE applies an extended Gaussian mixture model [17] to cluster SVGs, utilizing the same Gaussian‐process‐based prior employed previously during SVG detection. Another method SPARK‐X [18] identifies SVGs by testing the independence between two spot similarity matrices: one based on the expression levels of each gene and the other based on kernel‐transformed spatial locations. A gene is identified as an SVG if the null hypothesis of independence is rejected. Similar to SpatialDE and SPARK, SPARK‐X’s effectiveness relies on the ability of the chosen kernel to capture spatial patterns. Following SVG detection, both SPARK and SPARK‐X apply a log transformation to the raw count data, scaling the counts of each SVG to zero mean and unit variance across spots, and then perform hierarchical agglomerative clustering to group SVGs into co‐expression modules. MERINGUE [19], a graph‐based method, constructs a neighborhood graph of spatial spots using Delaunay triangulation, assigning binary edge weights (wij=1 if two spots are connected, wij=0 otherwise). It detects SVGs based on Moran’s I, transforming it into a z‐statistic and calculating a one‐sided p‐value for each gene. To cluster SVGs, MERINGUE calculates a spatial cross‐correlation index for each SVG pair, generating a spatial cross‐correlation matrix that is then used in hierarchical clustering to group SVGs into gene co‐expression modules. Although these methods incorporate clustering, their primary focus remains on identifying SVGs, and their clustering steps often rely on the same assumptions or data transformations used during SVG detection.
Certain spatial patterns may involve many co‐expressed genes, whereas others may involve only a few or even a single gene. Such pattern heterogeneity can result from factors such as pathological differences within tissue regions (e.g., tumor vs. normal areas in cancer) or the presence of diverse or rare cell types [16, 20]. Consequently, gene co‐expression modules may be imbalanced, with some containing many genes, whereas others are small or consist of a single gene. Standard clustering methods, such as mixture‐model‐based clustering [17] and hierarchical clustering [21], often lack the flexibility to capture such imbalanced clusters effectively. Although this challenge has been addressed in other contexts, such as time‐course gene expression [22, 23], it has not, to our knowledge, been specifically addressed in the context of ST data.
We present spatial transcriptomics iterative hierarchical clustering (stIHC), a novel two‐step method designed to identify gene co‐expression modules, including those with imbalanced sizes. In the first step, stIHC models the expression levels of SVGs across the spatial domain using a generalized penalized regression framework [24, 25], enabling efficient modeling of spatial variation over complex domains, including those with intricate geometries. This step is essential because gene expression data typically exhibit substantial noise while also following an intrinsically smooth spatial pattern. By modeling a gene’s expression levels at spatial spots as a smooth, continuous function rather than discrete values, it becomes possible to evaluate gene similarity across the entire tissue while accounting for spatial dependencies. This approach integrates relationships between adjacent spatial spots rather than treating a gene’s expression levels at different spots as independent observations. In the second step, stIHC clusters SVGs into co‐expression modules based on the estimated model coefficients. Recognizing that standard clustering methods often struggle with imbalanced gene clusters and fail to capture rare spatial expression patterns, stIHC employs the recently developed functional iterative hierarchical clustering (funIHC) [23]. This approach is specifically designed to handle imbalanced clusters, ensuring robust performance in identifying gene modules of any size that exhibit similar spatial expression patterns across a tissue slice.
We evaluate the effectiveness of stIHC across three simulated and four ST datasets generated using 10x Visium, 10x Xenium, and Spatial Transcriptomics technologies. Its performance is compared against the clustering approaches used by SpatialDE [15], SPARK [16], SPARK‐X [18], and MERINGUE [19] for grouping SVGs into co‐expression modules. The performance of stIHC is assessed using the adjusted rand index (ARI) and the Davies–Bouldin index (DBI) to determine its ability to identify coherent spatial gene co‐expression modules and capture rare or unique spatial patterns, which are often missed by the clustering approaches employed in the four SVG detection methods. Additionally, we examine whether the identified co‐expression modules are not only spatially coherent but also functionally relevant by analyzing their biological annotations and evaluating alignment with the known roles of the corresponding anatomical regions.
The contributions of this paper are as follows. First, we introduce stIHC, a novel clustering method for ST data that identifies spatial co‐expression gene modules, including those with unique spatial expression patterns often overlooked or merged into larger clusters by existing methods. Second, our method is data‐driven and parameter‐lean, requiring no user‐defined input parameters. Finally, by revealing gene co‐expression modules, our method provides valuable insights into the spatial organization of gene expression, establishing a powerful framework for exploring tissue structure and spatially coordinated biological processes.
The remainder of the paper is organized as follows. Section 2.1 introduces our proposed clustering method, stIHC. Section 2.2 presents simulation studies evaluating the performance of stIHC compared to existing approaches. Section 2.3 explores the functional annotation of gene clusters identified by stIHC. Section 2.4 demonstrates the application of stIHC to ST data from the mouse olfactory bulb. Finally, Section 3 concludes with a summary and discussion.
RESULTS
2
The stIHC methodology
2.1
ST data comprise gene expression levels measured at n distinct spatial locations, or “spots,” within a tissue slice. For each spot j=1,…,n, its 2D spatial location is denoted as sj=sj1,sj2⊤∈R2. The expression level of gene i=1,…,G at spot j is represented by yji∈R. Our method, stIHC, identifies gene co‐expression modules through a two‐step process: (1) modeling gene expression across a 2D tissue slice using a generalized penalized regression framework [24, 25], and (2) clustering the resulting basis coefficients using funIHC, a functional iterative hierarchical clustering algorithm [23].
Step 1: Modeling gene expression in a 2D tissue slice
2.1.1
Consider n locations s=s1,…,sn within a 2D tissue slice Ω⊆R2. At location sj, the observed expression level of gene i is denoted by yji∈R. We assume the response variable yji follows a distribution within the exponential family with mean μi. The exponential family includes many common distributions, such as Poisson, Binomial, Gamma, and Normal [26]. Thus, this framework can be used for modeling raw gene expression counts (e.g., Poisson) or normalized gene expression levels (e.g., Normal).
To begin, the slice of interest Ω is partitioned into smaller regions using Delaunay triangulation [27, 28], an efficient algorithm for handling nonregular geometries, such as the irregular tissue shape shown in the left panel of Figure 1. Locally supported polynomial functions are defined over these triangles, providing a set of basis functions. Let ϕ(s)=ϕ1(s),…,ϕK(s) represent an n×K matrix containing K piecewise linear basis functions evaluated at the locations s. Using these basis functions, we model μi via a generalized penalized regression framework [24, 25].
Let g be a continuously differentiable and strictly monotone canonical link function, such as g(μ)=log(μ) (Poisson) or g(μ)=μ (Normal), and let f represent a smooth spatial field over Ω. Then:
where ci are the coefficients of the basis function expansion approximating the behavior of gμi.
The coefficients ci are estimated by maximizing a penalized log‐likelihood functional:
where l(⋅) is the log‐likelihood, λ>0 is a smoothing parameter, and the Laplacian Δfi=∂2fi∂s12+∂2fi∂s22 measures the local curvature of the spatial field fi. Increased values of the smoothing parameter λ yield smoother estimates of f, whereas decreased values of λ give estimates that more closely fit the data. The Laplacian serves as the optimal penalization choice for mitigating the effects of noise present in the data [24]. Nevertheless, if the user has domain‐specific insights that require penalizing deviations from a more complex partial differential equation, such adjustments are possible; see Ref. [29] for further details. Using a penalized iterative reweighted least squares algorithm [25], we solve the optimization problem in Equation (2) to estimate the coefficients ci for each gene i for a fixed λ. The optimal smoothing parameter λ is estimated by minimizing the generalized cross‐validation criterion, as proposed by Craven and Wahba [30], across all genes, thereby determining a unified λ that ensures consistent penalization throughout the entire gene set. Let C represent the G×K matrix of coefficients corresponding to the optimal λ configuration, where C=c1,…,cG⊤ denotes the complete set of coefficients for all G genes.
Step 2: Clustering of genes based on spatial expression patterns
2.1.2
Using the resulting K basis coefficients for all genes, the distance metric di,j, quantifying the dissimilarity between two genes i and j, is computed as one minus the Spearman correlation ρi,j between ci and cj, that is, the ith and jth rows of C. The clustering procedure proceeds as follows:
- Define αmin and αmax as the minimum and maximum values of the Spearman correlation computed across all possible gene pairs {ρi,j}i,j=1G. Construct a grid spanning αmin,αmax, comprising U equally spaced values. For each αu, where u=1,…,U:
-
Cluster: Perform hierarchical clustering on the genes using the distance metric di,j with a threshold of 1−αu. Use the average linkage method to determine the dissimilarity between clusters. Denote the resulting number of clusters as P.
-
Merge: Let μp represent the center of cluster p (p=1,…,P). Treat each μp as a new gene and apply the same criterion as in step i to decide whether any cluster centers μp should be merged. If merging is identified, consolidate the respective clusters into a single one, resulting in R clusters, where R≤P.
-
Prune: For each cluster r=1,…,R, if the Spearman correlation ρi,r between μr and the ith gene within falls below αu, remove gene i from cluster r. All removed genes are allocated into single‐gene clusters. This results in a total of R+S clusters, where S is the number of removed genes.
-
Repeat steps ii and iii: Continue merging and pruning until the cluster assignment converges.
-
Repeat step ii: Repeat the merging step until all between‐cluster‐center Spearman correlations are less than αu.
-
For the rationale behind steps iv and v, please refer to Appendix A.
-
2Determine the optimal value αopt from α1,…,αU by maximizing the average silhouette value across genes. The silhouette value for the ith gene is defined as follows:
where ai is the average distance from the ith gene to other genes in the same cluster, and bi is the minimum average distance from the ith gene to genes in any other cluster. The silhouette value ranges from −1 to 1, with higher values indicating stronger similarity within clusters and greater dissimilarity between clusters.
Simulations
2.2
We conducted three simulations to evaluate the clustering performance of stIHC for ST data in comparison to the clustering approaches used in four SVG detection methods: SpatialDE [15], SPARK [16], SPARK‐X [18], and MERINGUE [19]. Implementation details are provided in Section 4.2. These simulations were based on the 10x Visium sagittal mouse brain slice dataset, available in the Seurat R package [31]. The first simulation evaluates performance under an ideal scenario with balanced equally sized clusters (Section 2.2.1). The second simulation tests performance with imbalanced cluster sizes (Section 2.2.2). The third simulation examines performance on sparse imbalanced clusters to mimic the real ST data analyzed in Section 2.4 (Section 2.2.3). To assess the stability of stIHC, we repeated the simulations 100 times and observed that the method consistently produced identical clustering results across all runs. This demonstrates that stIHC is a robust method. We also note that stIHC does not introduce variability due to random initialization or stochastic procedure, ensuring reproducible clustering on a particular dataset. Moreover, Appendix B presents the execution time in seconds for each method and simulation scenario. stIHC consistently demonstrates efficient computational performance across the three simulation settings.
Clustering performance with equally sized spatial co‐expression modules
2.2.1
This first simulation evaluates clustering performance when all gene modules have clearly distinguishable spatial expression patterns and an equal number of genes in each module. To represent distinct spatial patterns, we selected four unique gene expression patterns from the 10x Visium sagittal mouse brain slice dataset, each serving as the representative pattern for a different module. For a visual representation of these spatial expression patterns, refer to the top row of Figure 2.
Four genes from the 10x Visium sagittal mouse brain slice dataset representing the four clusters (top row). Example of one simulated gene from each cluster generated using scDesign3 (bottom row).
We simulated the spatial expression patterns of 100 genes from the four representative patterns (25 genes generated per pattern) using scDesign3 [32]. For illustration, an example of one simulated gene from each cluster is displayed in the bottom row of Figure 2. Full details of the data simulation process are provided in Section 4.1.
We applied five clustering methods: our proposed method, stIHC, along with the clustering approaches in SpatialDE, SPARK, SPARK‐X, and MERINGUE. We evaluated the performance of each method using the ARI [33], which ranges from −1 to 1, with 1 indicating complete agreement between the clusters and the ground‐truth gene modules. All methods successfully partitioned the 100 genes into their respective clusters, achieving a perfect clustering solution with an ARI of 1. These results indicate that the methods effectively clustered genes with similar spatial expression patterns in this balanced scenario with distinct clusters.
Clustering performance with imbalanced spatial co‐expression modules
2.2.2
The first simulation provided a baseline assessment of performance under ideal conditions with balanced gene modules. However, real‐world ST data may have gene modules of varying sizes. To simulate this more realistic scenario, we adjusted the dataset from the first simulation to create imbalanced gene modules. Specifically, we simulated 6 genes in the first module with expressions resembling Calb2, 2 genes in the second module resembling Ttr, 16 genes in the third module resembling Gpr88, and 25 genes in the fourth module resembling Cck (see Figure 2). Performance was assessed using the ARI and the DBI [34], which ranges from 0,∞, with lower values indicating better clustering quality. Results are summarized in Table 1.
When applied to this imbalanced dataset, stIHC was the only method to correctly identify all four modules along with the correct number of genes in each module. It achieved an ARI of 1 and the lowest DBI of 0.49. In contrast, SPARK, SPARK‐X, and MERINGUE identified the same three clusters (sizes: 6, 18, and 25 genes) and failed to detect the smallest cluster of just two genes, resulting in an ARI of 0.94 and a DBI of 0.54. SpatialDE performed the worst, identifying four clusters but completely misrepresenting the spatial patterns, resulting in cluster sizes of 16, 13, 8, and 12 genes, with a poor ARI of 0.11 and a high DBI of 4.54. These results highlight the limitations of SPARK, SPARK‐X, MERINGUE, and SpatialDE, which tend to merge smaller clusters into larger ones. In contrast, stIHC effectively captures and preserves rare spatial patterns.
We further evaluated the clustering methods by comparing the mean spatial expression patterns of the clusters identified by each method to the ground‐truth patterns used in the simulation. As shown in Figure 3, methods that failed to capture the ground‐truth gene modules produced distorted spatial expression patterns. The ground‐truth spatial expression patterns of gene modules are presented in the top row, followed by the mean patterns identified by stIHC (second row), SpatialDE (third row), and SPARK, SPARK‐X, and MERINGUE (bottom row).
The ground‐truth mean spatial expression patterns for each gene module in the imbalanced simulation (top row), followed by the mean spatial expression patterns identified by stIHC (second row), SpatialDE (third row), and SPARK, SPARK‐X, and MERINGUE (bottom row).
Specifically, stIHC accurately identified the correct clusters, maintaining clear, consistent, and distinct mean spatial expression patterns for each cluster. It was the only method to correctly capture the spatial pattern of Ttr. SpatialDE partially captured the patterns of Calb2 and Cck but failed to capture the patterns of Ttr and Gpr88. SPARK, SPARK‐X, and MERINGUE failed to capture the spatial pattern of Ttr but adequately captured the patterns of Calb2, Gpr88, and Cck. These findings demonstrate that SpatialDE, SPARK, SPARK‐X, and MERINGUE are limited in their ability to capture the spatial expression patterns of rare or small gene modules, whereas stIHC effectively identifies these patterns.
Clustering performance with sparse spatial resolution and imbalanced spatial co‐expression modules
2.2.3
This simulation evaluates clustering performance on a sparse version of the imbalanced dataset described in Section 2.2.2, where 260 spatial locations were randomly sampled. This sparse spatial resolution mimics the real dataset analyzed in Section 2.4 and is designed to test stIHC’s effectiveness under reduced spatial resolution. We maintained the same imbalanced module sizes as in Section 2.2.2: 6 genes in the first module reflecting the spatial expression pattern of Calb2, 2 genes in the second module reflecting the pattern of Ttr, 16 genes in the third module representing Gpr88, and 25 genes in the fourth module representing Cck. Performance was assessed using the ARI and the DBI, with results summarized in Table 2.
When applied to this sparse imbalanced dataset, stIHC identified five clusters with sizes 6, 1, 1, 16, and 25. It correctly recovered the number of genes in the Calb2, Cck, and Gpr88 modules but split the smallest module (Ttr) into two separate clusters, each containing a single gene. Despite this, stIHC achieved the highest ARI of 0.99 and the lowest DBI of 0.37 among all methods examined. It is worth noting that stIHC was the only method capable of preserving the spatial expression pattern of the smallest module, as evidenced by comparing its mean spatial expression patterns to the ground truth in Figure 4.
The ground truth mean spatial expression patterns for each cluster in the sparse imbalanced simulation (top row), followed by the mean spatial expression patterns identified by stIHC (second row), SpatialDE (third row), and SPARK, SPARK‐X, and MERINGUE (bottom row).
In contrast, SPARK, SPARK‐X, and MERINGUE identified only three clusters with sizes of 6, 18, and 25 genes, consistent with the results from Section 2.2.2. These methods failed to detect the smallest module of two genes, leading to a lower ARI of 0.94 and a higher DBI of 0.54. Although the mean spatial expression patterns for the identified clusters retained the patterns of Calb2, Cck, and Gpr88, the spatial pattern of the smallest module (Ttr) was entirely missed, as shown in Figure 4.
SpatialDE performed the worst, identifying only two clusters with sizes 31 and 18. This resulted in a significantly lower ARI of 0.05 and a much higher DBI of 4.45. The mean spatial expression patterns produced by SpatialDE were significantly distorted compared to the ground truth, highlighting its limitations in handling sparse spatial resolution and imbalanced spatial co‐expression modules (Figure 4).
stIHC identifies functional gene modules by clustering spatial expression patterns
2.3
To further validate stIHC, we applied it to three real ST datasets: 10x Visium mouse brain, 10x Xenium mouse brain, and 10x Visium human lung cancer (Sections 2.3.1–2.3.3, respectively). Our goal was to determine whether genes with similar spatial expression patterns also share biological functions. Specifically, we evaluated the biological coherence and functional distinctiveness of each detected co‐expression module. For the 10x Visium mouse brain and 10x Xenium mouse brain datasets, we assessed the spatial expression pattern similarity within clusters, alignment with known tissue structure, and gene ontology (GO) enrichment analysis. We also examined the consistency of results across the two platforms and tissue orientations. For the 10x Visium human lung cancer dataset, we evaluated the similarity of spatial expression patterns within clusters and performed GO enrichment analysis.
Each dataset consisted of 50 curated genes identified as SVGs. This selection ensures high confidence in the SVGs and removes potential biases from SVG detection methods, allowing for a focused evaluation of clustering performance. Although our method supports datasets of any size, we limited this analysis to curated SVGs for consistency. Details of SVG selection are provided in Section 4.1.
10x Visium mouse brain dataset
2.3.1
The 10x Visium mouse brain dataset was preprocessed to include 50 curated SVGs (Section 4.1). Using stIHC, we identified five co‐expression modules containing 2, 10, 16, 7, and 15 genes, respectively. Figure 5A shows the mean spatial expression pattern for each cluster, along with two representative genes per cluster. These spatially distinct patterns demonstrate that stIHC effectively captures the inherent spatial structure of the data.
stIHC applied to the 10x Visium mouse brain dataset. (A) Mean spatial expression patterns of the five clusters with two representative genes per cluster. (B) GO enrichment terms for the clusters. (C) Anatomical regions from the Allen Mouse Brain Atlas and Allen Reference Atlas (mouse.brain‐map.org and atlas.brain‐map.org) [14, 35].
To evaluate the biological relevance of these clusters, we performed GO enrichment analysis using the clusterProfiler R package [36]. Each cluster was enriched for unique biological processes (Figure 5B). Metrics such as p.adjust (adjusted p‐value), gene ratio (the proportion of genes associated with a GO term), and count (the number of associated genes) confirm the distinct functional identities of the clusters. Full gene lists are provided in Appendix C.
We compared the spatial expression patterns of these clusters with anatomical regions from the Allen Mouse Brain Atlas (Figure 5C). Each cluster corresponded to a specific brain region.
- Cluster 1 (2 genes): Thalamus, aligning with pathways related to sensory and motor signal relay [37].
- Cluster 2 (10 genes): Hippocampus, including the dentate gyrus and Ammon’s horn, with GO terms highlighting memory processes [38, 39].
- Cluster 3 (16 genes): Hypothalamus, associated with hormone regulation and production [40].
- Cluster 4 (7 genes): Thalamus and hippocampus, enriched for catabolic processes.
- Cluster 5 (15 genes): Cortical subplate, linked to migration and ion transport [41].
These findings confirm that the spatially distinct gene modules identified by stIHC align with the functional organization of the tissue.
10x Xenium mouse brain dataset
2.3.2
The 10x Xenium mouse brain dataset was similarly processed to include 50 SVGs (Section 4.1). stIHC identified three co‐expression modules containing 21, 17, and 12 genes, respectively. Figure 6A shows the mean spatial expression patterns of the clusters alongside two representative genes per cluster. The clear spatial distinctions reflect the effectiveness of stIHC in capturing spatial structures, even with different sequencing technologies and tissue orientations.
stIHC applied to the 10x Xenium mouse brain dataset. (A) Mean spatial expression patterns of the three clusters with two representative genes per cluster. (B) GO enrichment terms for the clusters. (C) Anatomical regions from the Allen Reference Atlas—Mouse Brain (atlas.brain‐map.org) [35].
GO enrichment analysis revealed distinct biological functions for each cluster (Figure 6B), aligning with anatomical regions from the Allen Mouse Brain Atlas (Figure 6C).
- Cluster 1 (21 genes): Hypothalamus, enriched for hormone regulation processes [40].
- Cluster 2 (17 genes): Hippocampus, associated with learning and cognition [38, 39].
- Cluster 3 (12 genes): Thalamus, linked to sleep regulation [42].
Despite differences in spatial resolution, molecular profiling, and tissue preparation between the 10x Visium and 10x Xenium platforms, stIHC consistently identified key regions (hippocampus, thalamus, and hypothalamus) with highly similar biological functions from both datasets. GO enrichment results for the hippocampus and hypothalamus clusters aligned across both platforms, reflecting learning and memory processes in the hippocampus and hormone regulation in the hypothalamus. These findings underscore the robustness and applicability of stIHC across ST platforms, demonstrating its ability to uncover biological phenomena beyond platform‐specific artifacts.
10x Visium human lung cancer dataset
2.3.3
The 10x Visium human lung cancer dataset underwent similar processing to include 50 SVGs (Section 4.1). Using stIHC, two co‐expression modules were identified, consisting of 16 and 34 genes, respectively. Figure 7A presents the mean spatial expression pattern for each cluster, accompanied by two representative genes from each cluster. These distinct spatial patterns illustrate that stIHC successfully captures the underlying structure of the data. We performed GO enrichment analysis and identified unique biological functions associated with each cluster (Figure 7B).
- Cluster 1 (16 genes): Associated with primary lung functions, as indicated by enrichment in muscle contraction, muscle system processes, and respiratory gaseous exchange by the respiratory system. Additionally, it is linked to the immune response, as evidenced by enrichment in the B cell receptor signaling pathway, which plays a role in the tumor microenvironment [43]. Furthermore, it includes genes involved in the cellular response to interleukin‐6 (IL‐6), which has been associated with lung cancer progression, resistance to antitumor therapies, and poor survival in lung cancer patients [44].
- Cluster 2 (34 genes): Associated with pathways that support tumor growth and development, including peptide hormone processing, which plays a known role in lung cancer development [45]; hemostasis, as dysregulation of the hemostatic system is common in lung cancer and serves as a prognostic indicator [46]; and tissue homeostasis, which contributes to tumor growth and progression [47].
stIHC applied to the 10x Visium lung cancer dataset. (A) Mean spatial expression patterns of the two clusters with two representative genes per cluster. (B) GO enrichment terms for the clusters.
The region in the top right corner, where Cluster 1 is highly expressed and Cluster 2 is lowly expressed, likely corresponds to a functionally distinct lung region enriched in primary lung functions and immune‐related processes. This aligns with the biological role of Cluster 1 genes in muscle contraction, respiratory function, and immune processes, suggesting that this area may represent non‐tumor tissue or a microenvironment supportive of normal lung function, in contrast to tumor‐associated regions where Cluster 2 is more highly expressed. These results indicate that the spatially distinct gene modules identified by stIHC are associated with unique biological functions. Furthermore, stIHC demonstrates the capability to uncover spatial gene co‐expression modules with potential disease‐specific relevance.
Application of stIHC to mouse olfactory bulb data
2.4
To demonstrate the practical application of stIHC, we analyzed ST data from replicate 11 of the mouse olfactory bulb, as described in Ref. [48]. This dataset, widely studied in previous research [15, 16, 19], is publicly available at The Spatial Research Lab’s website. Figure 8 shows a hematoxylin and eosin‐stained brightfield image of the mouse olfactory bulb, providing a structural reference for the tissue. In this analysis, we applied two preprocessing pipelines and SVG detection methods—MERINGUE and SpatialDE—to assess the flexibility and robustness of stIHC across different analytical workflows.
Hematoxylin and eosin‐stained brightfield image of the mouse olfactory bulb. Image from The Spatial Research Lab’s website (spatialresearch.org).
Data preprocessing and detection of SVGs
2.4.1
The raw dataset comprised 262 probe spots and 15,928 genes. Data preprocessing and SVG identification were conducted using the guidelines of MERINGUE and SpatialDE:
- MERINGUE: Probe spots with fewer than 100 detected genes were excluded, retaining 260 spots. Genes with fewer than 100 total reads were removed, resulting in 7365 genes. Raw counts were normalized to counts per million (CPM) values.
- SpatialDE: Genes with fewer than three total counts were filtered out. Variance stabilization was performed using Anscombe’s transformation to meet SpatialDE’s assumption of normally distributed data. Library size differences were adjusted by regressing out total counts for each gene. After preprocessing, 14,859 genes and 260 spots remained.
MERINGUE detected 886 SVGs, whereas SpatialDE identified 67 SVGs, highlighting substantial differences in SVG detection rates. Although this study focuses on these two methods, stIHC is compatible with any SVG detection method and adaptive to diverse pipelines.
Clustering SVGs into co‐expression modules
2.4.2
The log‐transformed counts (log(1+x)) of the detected SVGs were input into stIHC to cluster the genes into co‐expression modules based on their spatial expression patterns. Unlike methods such as SpatialDE, SPARK, and SPARK‐X, which require users to predefine the number of clusters, stIHC determines the optimal number of clusters automatically using an internal metric.
For the 886 SVGs detected by MERINGUE, stIHC identified three clusters containing 523, 341, and 22 genes. For the 67 SVGs detected by SpatialDE, stIHC identified two clusters containing 38 and 29 genes.
Figure 9 illustrates the spatial expression patterns of the three clusters identified by stIHC from the 886 SVGs detected using MERINGUE. The top row displays a representative gene from each cluster, whereas the middle and bottom rows illustrate the mean spatial expression patterns in the original spatial resolution and as a smoothed surface, respectively.
- Cluster 1: Genes in this cluster exhibited high expression localized in the granular cell layer, consistent with prior findings [16].
- Cluster 2: This cluster included genes with elevated expression in the central region but lacked a distinct spatial pattern.
- Cluster 3: Genes in this cluster showed high expression at the edges of the bulb, corresponding to the glomerular layer [16].
Three co‐expression modules identified by stIHC from 886 spatially variable genes (SVGs) detected by MERINGUE. Top row: representative genes for each cluster. Middle row: mean spatial expression patterns in the original resolution. Bottom row: smoothed spatial surfaces.
The spatial patterns of Clusters 1 and 3 align with findings reported in Refs. [15, 16, 19].
Figure 10 illustrates the two co‐expression modules identified by stIHC from the 67 SVGs detected by SpatialDE:
- Cluster 1: Genes in this cluster displayed high expression in the inner region, corresponding to the granular cell layer (analogous to MERINGUE’s Cluster 1).
- Cluster 2: This cluster showed high expression at the bulb’s edges, aligning with MERINGUE’s Cluster 3 (glomerular layer).
Two co‐expression modules identified by stIHC from 67 spatially variable genes (SVGs) detected by SpatialDE. Top row: representative genes for each cluster. Middle row: mean spatial expression patterns in the original resolution. Bottom row: smoothed spatial surfaces.
The absence of a cluster corresponding to MERINGUE’s Cluster 2 suggests that the larger SVG set identified by MERINGUE may include genes with less distinct or noisy spatial patterns, as Cluster 2 does not exhibit a clear spatial structure.
These results highlight the flexibility and robustness of stIHC in identifying meaningful spatial co‐expression modules from SVGs detected using different methods. Despite variations in preprocessing pipelines, the number of detected SVGs, and the resulting number of co‐expression modules, stIHC consistently uncovered biologically relevant spatial patterns. It is worth noting that stIHC preserved key spatial features, such as the granular and glomerular layers, across both MERINGUE and SpatialDE analyses, demonstrating its effectiveness in capturing the functional organization of the mouse olfactory bulb.
DISCUSSION AND CONCLUSION
3
We presented a novel clustering method, stIHC, designed to identify spatial co‐expression gene modules from ST data. stIHC demonstrated significant advantages over existing methods in simulations, especially in detecting rare and unique spatial co‐expression modules often overlooked by standard clustering approaches. Although traditional methods perform adequately in datasets with balanced cluster sizes, they falter when clusters are imbalanced. By accommodating imbalanced clusters, stIHC consistently outperformed other methods, successfully identifying spatial co‐expression modules even when represented by only a few genes.
When applied to ST datasets from the 10x Visium and 10x Xenium mouse brain, as well as the 10x Visium human lung cancer dataset, stIHC effectively identified biologically meaningful gene co‐expression modules. These modules exhibited distinct spatial expression patterns and, in the case of mouse brain datasets, corresponded to known anatomical regions of the brain, supporting the hypothesis that spatially organized gene expression reflects the tissue’s functional structure. GO enrichment analysis confirmed that genes within the same co‐expression module shared coherent and distinct biological functions, further reinforcing the biological relevance of spatially co‐expressed genes.
We evaluated stIHC’s robustness and generalizability across multiple datasets and ST technologies, including 10x Visium, 10x Xenium, and Spatial Transcriptomics. As demonstrated in Section 2.3, stIHC consistently identified spatially co‐expressed gene modules aligned with known biological pathways and anatomical structures, regardless of the technology or tissue orientation.
As shown in Section 2.4, stIHC can be integrated with any SVG detection method, offering a data‐driven and parameter‐lean clustering approach that eliminates the need for user‐defined input parameters. Its ability to preserve spatial patterns while detecting biologically and functionally relevant clusters makes it a powerful tool for exploring the complex spatial organization of tissues. The R implementation of stIHC is freely available on GitHub (CatherineH1/stIHC).
MATERIALS AND METHODS
4
Datasets
4.1
The 10x Visium sagittal mouse brain slice dataset is available in the Seurat R package. To represent distinct spatial patterns, we selected four distinct gene expression patterns, each serving as the representative pattern for a different cluster. Using the scDesign3 R package we simulated 25 copies of each gene, resulting in a total of 100 genes grouped into four clusters. We followed the tutorial available on GitHub (songdongyuan1994.github.io/scDesign3/docs/articles/scDesign‐spatial‐vignette) to generate this data.
The 10x Visium, 10x Xenium mouse brain data, and 10x Visium human lung cancer data used in Section 2.3 are available on GitHub (pinellolab/SVG_Benchmarking). scDesign3 was used to generate biologically realistic data with various spatial variability using real‐world ST data as a reference. The marginal distribution of expression for each gene was modeled using the function, fit_marginal (mu_formula = “s(spatial1, spatial2, bs = ‘gp’, k = 500)”, sigma_formula = “1”, family_use = “nb”), which fitted the data with a generalized GP model under the negative binomial distribution. The joint distribution of genes was modeled using the function, fit_copula (family_use = “nb”, copula = “gaussian”). Next, they extracted the mean parameters for each gene across all spots, denoted by μs(s) to remove spatial correlation between the spots, generating a nonspatially variable mean function μns(s), and used the function “sim_new” to generate simulation data as follows:
where α denotes the fraction of spatial variability in simulated gene expression. For our analysis in this paper we selected the genes for which α=1, meaning the expressions were generated from the GP model with the same spatial variability as the reference data.
The mouse olfactory bulb data for the replicate 11 sample can be downloaded from The Spatial Research Lab’s website (spatialresearch.org).
Existing method implementation
4.2
To cluster genes, SPARK and SPARK‐X use the hierarchical agglomerative clustering algorithm in the R package amap with the two optional parameters in the R function set to be Euclidean distance and Ward’s criterion, respectively. We used this function for a range of cluster values and selected the optimal clustering partition by maximizing the silhouette index.
MERINGUE was implemented by following the tutorial available at JEFworks‐Lab website (jef.works/MERINGUE/mOB_analysis).
The Bioconductor version of SpatialDE was implemented by following the tutorial available on Bioconductor (bioconductor.org/packages/release/bioc/vignettes/spatialDE/inst/doc/spatialDE.html).
AUTHOR CONTRIBUTIONS
Catherine Higgins: Formal analysis; funding acquisition; investigation; methodology; software; writing—original draft. Jingyi Jessica Li: Conceptualization; supervision; writing—review and editing. Michelle Carey: Funding acquisition; methodology; project administration; writing—review and editing.
CONFLICT OF INTEREST STATEMENT
The authors declare no conflicts of interest.
ETHICS STATEMENT
This article does not contain any studies with human or animal subjects performed by any of the authors.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Shah S , Lubeck E , Zhou W , Cai L . In situ transcription profiling of single cells reveals spatial organization of cells in the mouse hippocampus. Neuron. 2016;92(2):342–357.27764670 10.1016/j.neuron.2016.10.001PMC 5087994 · doi ↗ · pubmed ↗
- 210x Genomics . Inside visium spatial capture technology; 2019.
- 3Rodriques SG , Stickels RR , Goeva A , Martin CA , Murray E , Vanderburg CR , et al. Slide‐seq: a scalable technology for measuring genome‐wide expression at high spatial resolution. Sci Technol Humanit. 2019;363(6434):1463–1467.10.1126/science.aaw 1219 PMC 692720930923225 · doi ↗ · pubmed ↗
- 4Stickels RR , Murray E , Kumar P , Li J , Marshall JL , Di Bella DJ , et al. Highly sensitive spatial transcriptomics at near‐cellular resolution with Slide‐seqv 2. Nat Biotechnol. 2021;39(3):313–319.33288904 10.1038/s 41587-020-0739-1PMC 8606189 · doi ↗ · pubmed ↗
- 5Janesick A , Shelansky R , Gottscho AD , Wagner F , Williams SR , Rouault M , et al. High resolution mapping of the tumor microenvironment using integrated single‐cell, spatial and in situ analysis. Nat Commun. 2023;14(1):8353.38114474 10.1038/s 41467-023-43458-x PMC 10730913 · doi ↗ · pubmed ↗
- 6Chen KH , Boettiger AN , Moffitt JR , Wang S , Zhuang X . Spatially resolved, highly multiplexed rna profiling in single cells. Sci Technol Humanit. 2015;348(6233):aaa 6090.10.1126/science.aaa 6090 PMC 466268125858977 · doi ↗ · pubmed ↗
- 7Lubeck E , Coskun AF , Zhiyentayev T , Ahmad M , Cai L . Single‐cell in situ RNA profiling by sequential hybridization. Nat Methods. 2014;11(4):360–361.24681720 10.1038/nmeth.2892 PMC 4085791 · doi ↗ · pubmed ↗
- 8Yan G , Hua SH , Li JJ . Categorization of 34 computational methods to detect spatially variable genes from spatially resolved transcriptomics data. Nat Commun. 2025;16(1):1141.39880807 10.1038/s 41467-025-56080-w PMC 11779979 · doi ↗ · pubmed ↗