Ultra-precision deconvolution of spatial transcriptomics decodes immune heterogeneity and fate-defining programs in tissues

Domain: biology
Published: 2026-03-20
Editorial outcome: major_revision

Abstract
Elucidating the spatial organization and functional specialization of immune cells within complex tissues remains challenging. We present UCASpatial, an ultra-precision spatial transcriptomics deconvolution algorithm utilizing entropy-based weighting to accurately map cell subpopulations. Benchmarking confirms its superiority in identifying low-abundant cell subpopulations and distinguishing transcriptionally heterogeneous cell subpopulations. Applying UCASpatial to human colorectal cancer, we reveal that chromosome 20q gain in individual cancer clones orchestrates a T cell-excluded microenvironment, associated with HERV-H silencing and impaired type I interferon responses. In murine wound healing models, we reveal spatiotemporal dynamics distinguishing scarring from regenerative phenotypes. Specifically, we identify a pro-fibrotic community comprising Igfbp5+ chondrocytes, Cd36+ Gpnmb+ Il1b- macrophages, and Fmod+ fibroblasts in scarring-healing mice (C57BL/6). We further demonstrate that IL11-IL11RA signaling within this triad drives the pro-fibrotic community formation and limits regeneration.

Together, UCASpatial serves as a versatile tool for deciphering fine-grained cellular landscapes and exploring intercellular mechanisms in complex and dynamic microenvironments.

Introduction
The ability to elucidate the spatial organization of cells within the microenvironment is crucial for understanding cellular functions and exploring intercellular mechanisms that govern development, physiology, and pathogenesis 1, 2. Emerging spatial transcriptomics (ST) technologies offer promising avenues for achieving this; however, current spot-based techniques face a significant limitation in lacking single-cell resolution. Methods like Visium 3, Tomo-Seq 4, Slide-seq 5, and HDST 6 amalgamate RNA-seq measurements from multiple cells or cell fractions at each spot. As a result, these techniques struggle with the precise identification of individual cell types within the microenvironment. Consequently, annotating cell types at each spot becomes imperative for accurately capturing spatial cell organization 7, 8.

To address this challenge, various computational approaches have been proposed to deconvolute spatial transcriptomic profiles into cellular compositions based on gene expression profiles of each cell type identified through single-cell RNA sequencing (scRNA-seq), such as SPOTlight 9, RCTD 10, cell2location 11, CARD 12, stereoscope 13, and Spotiphy 14. These methods effectively delineate significant differences in cell type compositions across diverse anatomical or pathological layers in well-structured tissues such as the brain, heart, and lymphatic structures, even at fine-grained resolutions 11, 12, 15, 16.

In recent years, a growing body of research has emphasized the importance of characterizing the complex microenvironment, particularly the influx of immune cells that orchestrate dynamic spatial patterns, varying between tissue homeostasis and non-homeostasis conditions. This complexity is evident in scenarios such as tissue damage 17, 18 and the tumor microenvironment 19-23. However, several immune cells, particularly those infiltrating into tumors, are known for their sparse presence and low transcriptional activity 24, 25. Their discontinuous spatial distributions in these microenvironments limit the effectiveness of algorithms relying on spot-neighbor information. Additionally, the low abundance and the context-dependent phenotypic variation of heterogeneous immune cells, particularly T cells, pose challenges for conventional deconvolution algorithms 26. Furthermore, their scarcity makes them difficult to capture accurately using histological images, thus constraining the applicability of image-based algorithms 27, 28. Consequently, the precise characterization of complex and dynamic microenvironments remains a significant challenge.

To address the aforementioned challenge, we developed UCASpatial, an ultra-precision cellular deconvolution algorithm for spatial transcriptomics based on entropy-weighted regression. A unique advancement of UCASpatial is its formulation of the Shannon entropy-based weights to quantify the ability of each gene for indicating cell identity and incorporate this information in the deconvolution process. This approach improves the robustness and accuracy in identifying low-abundant cell subpopulations and distinguishing transcriptionally heterogeneous cell subpopulations in ST, thus allowing for more precise and unbiased deconvolution of the cellular landscape.

By applying UCASpatial to human colorectal cancer (CRC), UCASpatial has delineated the spatial patterns of immune cells within the tumor immune microenvironment (TIME). By integrating the fine-grained TIME landscape with spatially inferred copy number variations (CNVs), we characterize two functional multicellular communities associated with T cell-rich and T cell-excluded TIMEs at a clonal resolution. Further, we demonstrate that the chr20q-gain in CRC tumor cells is associated with reduced expression of the tumor-intrinsic human endogenous retrovirus subfamily H (HERV-H) and impaired type I interferon response, which orchestrates the T cell-excluded TIME phenotype. In the ear wound healing models of both regenerative MRL/MpJ (MRL) and scarring-healing C57BL/6 mice (B6 mice), we illustrate the spatiotemporal dynamics of immune cell subpopulations involved in the healing process and unveil the spatial segregation of macrophage subpopulations embedded within distinct cellular communities at the wound sites. Moreover, we identify a previously uncharacterized Igfbp5+ chondrocyte subpopulation that coordinates with Cd36+ Gpnmb+ Il1b- macrophage and Fmod+ fibroblast at the wound bed, establishing a three-cell-type community with pro-fibrotic functions linked to scarring in B6 mice. By blocking the Fmod+ fibroblast differentiation with the PPARγ agonist rosiglitazone, or by inhibiting the IL11-IL11RA axis within the community in B6 mice, we observed a reduction in the formation of the pro-fibrotic community and an improvement in the regenerative ability of ear wounds in B6 mice.

Results
UCASpatial: an entropy-weighted regression model for spatial deconvolution of cellular landscapes

UCASpatial deconvolutes ST data into cellular compositions by leveraging cell identity-specific gene expression profiles (cGEPs) derived from paired or unpaired scRNA-seq data (Fig. 1). First, UCASpatial extracts cGEPs from the referenced scRNA-seq dataset with a set of user-provided cell-subpopulation labels (fig. S1, A to C). It strategically selects cells meeting a meta-purity threshold to enhance discrimination of referenced cell identities, and then extracts cGEPs through optimized Non-negative Matrix Factorization (NMF). Second, UCASpatial formulates Shannon entropy-based weights for each gene within the extracted cGEPs to enhance genes indicative of cell identities (Methods). Third, UCASpatial employs Weighted Non-Negative Least Squares (WNNLS) to deconvolute the mRNA counts of each spot in the ST data using the weighted cGEPs, thereby estimating the cellular compositions. Additionally, UCASpatial offers a suite of downstream analysis tools. These include identifying tumor clones and their corresponding tumor immune microenvironment (TIME), as well as spatially resolved quantification of transposable elements (TEs). To ensure versatility in applications, UCASpatial is seamlessly integrated with the Seurat framework.

Benchmarking the performance of UCASpatial in the simulated ST dataset

To assess UCASpatial’s performance, we generated a simulated ST dataset comprising 750 ST spots designed to mimic the real cellular composition of the tumor microenvironment. Each spot was constituted by the transcripts of 5 to 20 cells sampled from scRNA-seq data of human CRC, encompassing 20 reference cell subpopulations (fig. S1, A and B, and Methods). The cellular composition of simulated spots was designed to recapitulate the spatially correlated cell-type distribution patterns observed in tumor microenvironments, particularly the relationship between epithelial, stromal, and immune cell frequencies across tumor core, margin, and stromal regions (fig. S1D and Supplementary Table 1) 29, 30.

To assess the performance of UCASpatial on the simulated ST dataset with varying complexity, we ranked the complexity of these simulated spots based on their cellular composition entropy. Subsequently, we categorized the simulated spots into three groups (low, medium, and high complexity), where higher entropy denotes increased complexity (Fig. 2A and fig. S1D). We compared UCASpatial with several recently established methods, including SPOTlight, RCTD, cell2location, CARD, stereoscope, and Spotiphy, using Root Mean Square Error (RMSE) to first quantify the deviation between estimations and ground truth at the spot level (Methods).

UCASpatial exhibited significantly superior performance across all complexity levels (median RMSE = 0.051) compared to the aforementioned methods, with 10%, 21%, 31%, 31%, 36%, and 35% improvement in terms of RMSE compared to RCTD (0.058), SPOTlight (0.066), CARD (0.075), Spotiphy (0.075), cell2location (0.081), and stereoscope (0.080), respectively (Fig. 2B and fig. S1E). Additionally, UCASpatial demonstrated a comparable Pearson Correlation Coefficient (PCC) to cell2location and significantly outperformed other methods, when comparing estimations with the ground truth (fig. S1F). Moreover, we also assessed the ability for identifying the existence of cell subpopulations within spots using the F1 score. Consistently, UCASpatial significantly outperformed other methods across all complexity levels with the highest F1 score (with a 6-25% improvement) (Fig. 2B and fig. S1G).

Additionally, to assess the robustness of UCASpatial across different cell numbers, we categorized spots into three groups based on the number of cells per spot: 5-9 cells, 10-12 cells, and 13-20 cells (fig. S1H). Our analysis revealed that UCASpatial maintained robust performance across these varying cell numbers, whereas the performance of cell2location and stereoscope dropped as cell number increased (fig. S1I). These results highlight UCASpatial’s robustness in handling diverse cell densities, thereby confirming its adaptability and reliability.

We next evaluated the precision and robustness of UCASpatial in deconvoluting different cell subpopulations. Through in-silico analyses using three distinct reference scRNA-seq labels with varying numbers of subpopulations (7, 11, and 20), we observed that the advantages of UCASpatial became more pronounced as the number of subpopulations increased compared with other methods (fig. S1, J and K). With 20 subpopulations, UCASpatial exhibited superior performance across all metrics, including RMSE, PCC, and F1 score, highlighting its proficiency in deciphering fine-grained spatial landscapes (fig. S1K). Specifically, examination of PCC between estimations and ground truth for each method under 20 subpopulations condition unveiled UCASpatial’s superior performance, achieving the highest PCC in 15 out of 20 cell subpopulations (Fig. 2C), particularly notable for NK cell and CD8+ T cell subpopulations (with a 39-187% improvement) (Fig. 2D and fig. S1L). By categorizing the PCC values of each cell subpopulation into four groups, we noticed UCASpatial exhibited PCC values surpassing the medium-high correlation level (PCC > 0.5) for 19 out of 20 cell subpopulations, with the sole exception CD4+ Tfh (PCC = 0.47) while remaining the best compared to alternative methods (Fig. 2, C and E). In contrast, other methods yielded PCC values below the medium-high correlation level for at least 30% of cell subpopulations, with some even falling below 0.25 (Fig. 2E and fig. S1M). This comparison underscored UCASpatial’s capability in discriminating cell subpopulations, even with subtle transcriptomic differences (fig. S1N). Furthermore, assessments using the PCC and F1 score reaffirmed UCASpatial’s precise and robust decomposition across each cell subpopulation (fig. S1, O and P).

Considering the challenge of ST analysis in precisely mapping the spatial location of low-abundant cell subpopulations, we next investigated UCASpatial’s sensitivity in annotating cell subpopulations with low frequencies within spots. We categorized the proportion of each cell subpopulation within simulated spots into ten groups, forming a gradient that spanned from 0-10% to 100% (Fig. 2F). Subsequently, we calculated the RMSE across various proportions to evaluate the deviation of cell-subpopulation proportions. We observed that UCASpatial exhibited the lowest overall RMSE (averaged by each cell subpopulation, medium RMSE = 0.051) across different proportions with a 7-31% improvement over the other methods (Fig. 2G). Moreover, we specifically examined the RMSE of different cell subpopulations, revealing that UCASpatial notably outperformed other methods for cell subpopulations with subtle transcriptomic differences such as CD8+ Tex and CD8+ Tem (with 25-63% and 13-26% improvements, respectively) (Fig. 2G and fig. S1Q). Consistent results were observed in assessments using the PCC and F1 score (fig. S1R), highlighting UCASpatial’s superiority in accurately annotating cell subpopulations with low frequencies within spots, even with subtle transcriptomic differences.

Furthermore, we evaluated the performance enhancement afforded by UCASpatial’s entropy-based weighting and meta-purity filtering designs (fig. S1, S to U). We confirmed that the meta-purity filtering process did not result in the loss of any predefined cell subpopulations, with PCC values before and after filtering consistently exceeding 0.9 (fig. S1, V to X). We found that the entropy-based weighting design showed improvements of 22%, 50%, and 17% in terms of RMSE, PCC, and F1 score, respectively, compared to UCASpatial without weighting. The meta-purity filtering design exhibited improvements across different parameters with 4-11%, 1-8%, and 1-3% in terms of RMSE, PCC, and F1 score. In addition, the same analyses as before showed that UCASpatial, without entropy-weighting and meta-purity filtering, lost the aforementioned improvements (fig. S1, E to R).

We evaluated the runtime and memory usage of UCASpatial and other methods (Methods), analyzing the impact of the number of spots in spatial data and the number of cells in scRNA-seq data on computational resources. UCASpatial showed stable time costs, processing data in under 1 hour as the number of spots increased (fig. S1Y). Cell2location encountered memory errors on our GPU platform with the largest dataset (50,000 spots and ~12,500 cells). Other methods, except stereoscope, were minimally affected by the scRNA-seq reference size (fig. S1Y). Using a simulated dataset (50,000 spots, ~12,500 cells, 20 cell subpopulations), UCASpatial demonstrated the lowest memory usage (< 26GB) and highest efficiency (fig. S1Z). Collectively, our in-silico benchmarking demonstrates that UCASpatial enables precise and sensitive spatial deconvolution, excelling in discerning cell subpopulations with subtle transcriptomic differences, even in low proportions.

UCASpatial demonstrates broad applicability and high precision across multi-resolution platforms

To test whether UCASpatial is applicable across platforms with varying resolutions, we applied UCASpatial to 10x Visium HD data from a colorectal cancer (CRC) sample 31, using the recommended bin size of 8 μm × 8 μm (fig. S2A and Methods). We compared UCASpatial with RCTD, a method also suitable for processing 10x Visium HD data using PCC between estimated proportions of cell subpopulations and marker-gene scores in the absence of ground truth 16. We observed that UCASpatial achieved the highest PCC, and exhibited over 20 times faster than RCTD on the runtime evaluation (fig. S2, B to D, Methods).

We also analyzed Slide-seq V2 data from the mouse hippocampus to evaluate UCASpatial’s performance at single-cell resolution (fig. S2, E to H). Since strictly defined ground-truth cell type labels are not available experimentally, we adopted the widely accepted evaluation strategy established by the CARD paper 12. This involves verifying whether the deconvolved cell types accurately reconstruct known anatomical structures (e.g., CA1, CA3, Dentate Gyrus). Specifically, we visualized the predicted top 1 cell type for each location (fig. S2F) and quantitatively evaluated the concordance between method-predicted cell-type proportions and known region-specific marker expression in three anatomical structures 12 by calculating PCC (fig. S2G). Our results show that UCASpatial is highly concordant with other methods, successfully captures key anatomical structures, and confirms a comparable level of accuracy and spatial resolution on this neurobiological dataset (fig. S2, G and H).

Together, these results demonstrate UCASpatial’s broad applicability and high precision across platforms with varying resolutions, from bin resolution (10x Visium HD) to single-cell resolution (Slide-seq V2).

Benchmarking the performance of UCASpatial using Visium HD-derived pseudo-spots

To obtain synthetic data that better capture spatial correlation structures between neighboring spots and cell-subpopulation-specific capture efficiency variations, we generated synthetic spots by spatially aggregating segmented single-cell transcriptomes from three 10x Genomics Visium HD colorectal cancer tissue sections31 into 49 µm × 49 µm bins (fig. S2, I to K), closely matching standard Visium spot size. This approach preserves local spatial correlations and yields 39,813 pseudo-spots with ground-truth cell-subpopulation compositions reflecting the tumor microenvironment’s complexity.

Benchmarking across this dataset recapitulated prior performance trends. UCASpatial achieved the highest average PCC (0.87) and F1 score (0.82), alongside the lowest RMSE (0.05) among other methods at the spot level (fig. S2L), even under scenarios of varying complexity (fig. S2M). Cell-subpopulation-level evaluation further confirmed UCASpatial’s superior accuracy in deciphering different cell subpopulations (fig. S2, L and N). This data further supports UCASpatial’s accurate and stable performance on biologically realistic synthetic data.

UCASpatial accurately annotates cell subpopulations in human CRC TIMEs

Spatial omics has provided a powerful tool for dissecting the intricate interplays between tumor alterations and the co-evolving immune microenvironment 21, 32, 33. However, the characterization of immune cell compositions, particularly those infiltrating into distinct tumor clone regions, remains a formidable challenge due to their low abundance and limited transcriptional output 21, 25, 34. To investigate whether UCASpatial can precisely resolve the tumor immune microenvironment (TIME), we curated an ST dataset from eight tumor slices of seven patients with human CRC using the 10x Genomics Visium platform as the discovery cohort 35 (Fig. 3A, fig. S3A, and Supplementary Table 2). We leveraged an integrated scRNA-seq dataset comprising an in-house scRNA-seq with three patients, and a public scRNA-seq dataset with five patients (P1-5) 35, encompassing 34 cell subpopulations as a reference (Fig. 3, A and B, fig. S3, B to D, and Methods).

We annotated the spatial distribution of 34 cell subpopulations in human CRC ST discovery cohort using UCASpatial. Notably, stromal cell subpopulations showed consistent distribution with anatomical locations as expected from prior literature 36, 37, such as smooth muscle cells of the muscular layer and endothelial cells around blood vessels, along with pericytes (fig. S3E). Similarly, epithelial cell subpopulations displayed enrichment in the epithelial region, consistent with corresponding H&E staining and pathological annotations (fig. S3F). Upon further examination of immune subpopulations, we noticed B cell and various T cell and dendritic cell (DC) subpopulations, including CD4+ Tn, CD4+ Tfh, CD8+ Tex, CD8+ Tem, cDC2, and CCR7+ DC, formed a multicellular community resembling a tertiary lymphoid structure (TLS) (fig. S3G, see Supplementary Table 3 for abbreviation information) 38, 39. Additionally, we noticed the spatial exclusivity of HES1hi C1QChi macrophages and TREM2+ PDPN+ macrophages (fig. S3H), indicating their distinct functions in TIMEs.

For comparative analysis, we also annotated this ST dataset using other methods. Remarkably, UCASpatial achieved the highest average PCC with a significant improvement of 17-104% over other methods between the cell-subpopulation proportions and signature expressions, especially in T cell subpopulations (with at least a 26% improvement) (fig. S3, I and J, Supplementary Data 1, and Methods). In alignment with in-silico benchmarking results, all evaluated methods correctly annotated structural cell subpopulations, for example, FAP+ myofibroblast (fig. S3J), yet UCASpatial showed a higher sensitivity in annotating immune cell subpopulations (fig. S3, I and J). For instance, UCASpatial performed better than other methods in identifying the spatial distribution of T cell subpopulations in almost all five TLS regions (fig. S3K). While UCASpatial depicted 4-6 T cell subpopulations consistent with their marker gene expressions, all compared methods only identified 0-2 T cell subpopulations (fig. S3L).

To assess the accuracy of immune cell annotation in the real dataset, we constructed the validation cohort. We first performed ST on an additional human CRC tumor sample (CRC8) and validated the ST results using mIHC staining on consecutive slices from the same tissue block utilized for ST (Fig. 3A, fig. S4, A to C, and Methods). The ST data were annotated using unpaired scRNA-seq data obtained from independent CRC datasets (Fig. 3B) 35. Quantitative analysis of cellular proportions via mIHC revealed that among 12 validated immune subpopulations, the majority of them had an average proportion of less than 10% within the pseudo-spots (Methods), with T cell subpopulations exhibiting the lowest proportions (fig. S4D). In alignment with the results from mIHC, UCASpatial correctly annotated the spatial distribution of these cell subpopulations (Fig. 3C and fig. S4, E to H). We further compared the cellular proportions estimated by other computational algorithms with those quantified by mIHC in tumor regions. Our analysis demonstrated that UCASpatial consistently achieved the highest PCC and the lowest RMSE for these immune cell subpopulations, particularly for T cell subpopulations, in tumor regions (Fig. 3C and fig. S4, E to HA). These findings underscore the capability of UCASpatial in accurately identifying low-abundant cell subpopulations, even in complex tissue contexts.

We found that the spatial distribution patterns of three CD8+ T cell subpopulations at the tumor-invasive front were distinct: CD8+ Tex and CD8+ Trm were more proximal to the tumor, whereas CD8+ Tem was relatively distal (Fig. 3, D and E). Also, there was a higher proportion of CD8+ Tex and CD8+ Trm infiltrated into the tumor region than CD8+ Tem (Fig. 3E). To validate these findings, we performed mIHC on human CRC samples from two additional patients (CRC9 and CRC10) in the validation cohort (Fig. 3A) and quantified the distance of each CD8+ T cell subpopulation to the tumor. This analysis revealed that PD-1+ CD8+ T cells and CD103+ CD8+ T cells were positioned closer to the tumor compared to GZMK+ CD8+ T cells (Fig. 3, F and G). Further calculation of the proportion of each CD8+ T cell subpopulation within the tumor region supported that PD-1+ CD8+ T cells and CD103+ CD8+ T cells exhibited higher tumor infiltration proportions compared to GZMK+ CD8+ T cells (Fig. 3, F and G). Additionally, we noticed UCASpatial without entropy-weighting and meta-purity filtering was unable to either identify T cell subpopulations in TLS regions or discriminate the spatial distribution of tumor-infiltrating CD8+ T cell subpopulations (Fig. 3D, fig. S3, I to L, and fig. S4, E to H). These results demonstrate that UCASpatial can accurately resolve cell subpopulations in the complex TIME.

Identification of functional multicellular communities associated with T cell-rich and T cell-excluded clonalTIMEs in CRC

Tumor-intrinsic genomic alterations have collateral effects that modify the spatial organization of the TIMEs 40-42. To spatially characterize diverse TIMEs at a detailed clonal resolution, we generated an integrated immunogenomics platform termed ‘UCASpatial-clonalTIME’ to spatially link tumor genomics with the TIMEs. This platform enables the identification of tumor clones within tissue slices by inferring spatial copy number variations (CNVs) 33 and concurrently estimating the composition of the TIME for each clone (clonalTIME) (Fig. 3H, fig. S4I, and Methods). By assessing the proportions of CD4+ and CD8+ T cell subpopulations, we categorized the clonalTIMEs into two phenotypes: T cell-rich and T cell-excluded (Fig. 3I, fig. S4J, and Methods). We detected three out of thirteen tumor clones displaying the T cell-rich clonalTIME phenotype (Fig. 3I and fig. S4, J to K). Within these T cell-rich clonalTIMEs, all three subpopulations of CD8+ T cells—CD8+ Tem, CD8+ Trm, and CD8+ Tex—were notably enriched (Fig. 3I). Analysis of the CD4+ T cell compartment showed significant enrichments of CD4+ Tm and Th1/Th17 within these T cell-rich clonalTIMEs (Fig. 3I). Conversely, within T cell-excluded clonalTIMEs, we found a depletion of all T cell subpopulations. Instead, there was an enrichment of classical monocytes (cMo), neutrophils, and TREM2+ PDPN+ macrophages (Fig. 3I).

Immune cells within the TIMEs engage in crosstalk with non-immune cells that coordinately form multicellular communities. To identify these communities associated with two clonalTIME phenotypes, we calculated the PCC between the neighborhood-averaged proportions of various cell subpopulations including 22 immune subpopulations and 7 non-immune subpopulations across all spots from the 13 clonalTIMEs (Methods). Through hierarchical clustering followed by functional enrichment analysis (Methods), we observed three multicellular communities discriminated by different cellular compositions and specialized functionalities: 1) Community1 encompassed HES1hi C1QChi macrophages, CD8+ Tem, CD8+ Tex and inflammatory CD34+ fibroblast (CD34+ iFib) (Fig. 3J and fig. S4M), with enriched pathways in leukocyte mediated cytotoxicity and homeostasis (fig. S4N); 2) Community2 comprised endothelial cells, pericytes, TREM2+ PDPN+ macrophages, TREM2+ APOC1hi macrophages, FAP+ myofibroblasts, smooth muscle cells, and glial cells (Fig. 3J), with functional pathways related to angiogenesis, wound healing, and WNT signaling (fig. S4N); 3) Community3, including neutrophils and cMos (Fig. 3J), exhibited enriched pathways in cytokine-mediated signaling pathway and acute inflammatory response (fig. S4N).

To link these communities with T cell-rich versus T cell-excluded clonalTIME phenotype, we examined the score of each community in the respective clonalTIME phenotypes (Methods). We observed a significantly higher score of community1 within T cell-rich clonalTIMEs, while the score of community2 was higher in T cell-excluded clonalTIMEs (Fig. 3K). The score of community3 showed no significant differences between the two clonalTIME phenotypes (Fig. 3K). Moreover, an opposite enrichment trend of community1 and community2 was detected across 13 clonalTIMEs (PCC = -0.51) (Fig. 3L). We next closely examined community distributions in a particular ST slice (CRC4) containing both T cell-rich (clone 1 and 2) and T cell-excluded (clone 3) tumor clones. The result revealed a loss of community1 coincident with an increase of community2 and community3 as the clonalTIME transitioned from a T cell-rich to a T cell-excluded phenotype (Fig. 3, M and N). Together, we spatially characterized the multicellular communities within T cell-rich and T cell-excluded TIME at a clonal resolution and identified two of these functional communities associated with T cell-rich versus excluded clonalTIME phenotypes in CRC.

Identification of tumor genetic alterations associated with clonalTIME transition in CRC

To investigate whether genetic alterations are associated with immune phenotypic transitions in tumor clones, we conducted a comparison of CNVs across clones (fig. S5A and Supplementary Table 4). Our findings revealed a distinct pattern among T cell-excluded clones, exhibiting unique CNVs that are absent in T cell-rich clones (Fig. 4A). Notably, we observed a frequent copy number gain on chromosomes 20q, 20p, and 7p (chr20q-gain, chr20p-gain, and chr7p-gain) within T cell-excluded clones (Fig. 4A). This observation prompted us to hypothesize that these chromosomal gains may arise during clonal evolution and promote the transition within the tumor microenvironment from a T cell-rich to a T cell-excluded state.

To further investigate the relationship between genetic alterations and the immunophenotype of tumor clones, we focused on the aforementioned CRC4 ST slice containing three clones exhibiting sequential clonal evolution. Clone1, characterized by chr1p-loss-2, served as a common ancestor (Fig. 4, B and C). It gave rise to two distinct progenies: clone2 that acquired chr8p-gain and chr8q-gain, and clone3 with chr20q-gain, chr20p-gain, chr7p-gain, and two additional CNVs (Fig. 4, B and C). Notably, while the evolution from clone1 to clone2 was accompanied by maintenance of a T cell-rich TIME (Fig. 4C), the evolution to clone3 was marked by a transition to the T cell-excluded phenotype, evidenced by a significant reduction in CD8+ Tex, CD8+ Tem, and HES1hi C1QChi macrophages infiltration within their clonalTIMEs (Fig. 3, I and N, and Fig. 4C). This transition suggests that the genetic changes accumulated in clone3, particularly the chr20q-gain and associated CNVs, may be critical in shaping a T cell-excluded TIME.

For further validation of the relationship between the genetic changes specific to CRC4-clone3 and the T cell-exclusive immunophenotype, we conducted an integrative analysis utilizing the TCGA-COAD dataset, which encompasses DNA copy-number variation (Affymetrix 6.0 microarrays) and bulk RNA sequencing data with 471 samples. We identified the existence of CRC4-identified CNVs and estimated the infiltration levels of various immune cell subpopulations based on their signature expressions derived from bulk RNA-seq for each patient. We then calculated the PCC between the presence of these CNVs and estimated cell subpopulations (fig. S5B), revealing significant correlations between the presence of clone3-specific CNVs and the reduction in the signature expression of various T cell subpopulations (Fig. 4D).

Notably, all of these five clone3-specific CNVs (chr20q-gain, chr20p-gain, chr7p-gain, chr7q-gain, and chr17p-loss) exhibited a concurrent occurrence pattern within the TCGA-COAD dataset (Fig. 4E). To determine the primary CNV responsible for the diminished T-cell infiltration, we isolated the individual impacts of each CNV in the TCGA-COAD dataset based on their presence (Methods). Our analysis revealed that only the group harboring chr20q-gain, prevalent in 57.9% of CRC patients (Fig. 4D), exhibited a significant reduction in the expression levels of T cell marker genes (CD3G, CD8A, and CD4) when compared to the group without chr20q-gain (Fig. 4F and fig. S5C). Since chromosomal instability (CIN) has been previously reported to be associated with decreased immune cell infiltration 43, we stratified the samples by aneuploidy levels to control for its potential effects (Methods). Patients with chr20q-gain showed significantly lower T cell marker expression than chr20q-WT patients in both low and intermediate aneuploidy groups (fig. S5D), demonstrating a CIN-independent role of chr20q-gain in driving T cell exclusion. These findings suggest that the presence of chr20q-gain may play a pivotal role in driving the transition of the clonalTIME phenotype from T cell-rich to T cell-excluded.

Tumor-intrinsic endogenous retrovirus silencing correlates with chr20q-gain related T cell-excluded clonalTIME formation

To investigate the potential mechanism of the immune phenotypic transitions in tumor clones, we performed differential gene expression analysis in clone3, which had acquired a chr20q gain, compared to clone1 and clone2 of CRC4. Through functional enrichment analysis, we identified enrichment for pathways related to translation machinery among the upregulated genes in clone3 (Fig. 4G). Conversely, pathways related to type I interferon (IFN) production were notably enriched among the downregulated genes, indicating a potential suppression of immune response (Fig. 4G and fig. S5, E to G). Notably, chemokine-coding genes responsible for cytotoxic CD8+ T cell recruitment, such as CXCL9 and CXCL10, and genes associated with antigen presentation, including TAP1, TAP2, HLA-A, HLA-B, and B2M 44, 45, exhibited reduced expression (Fig. 4G and fig. S5, E to G). This observation potentially elucidates the T cell-excluded phenotype of clone3.

It is well-established that the activation of the type I IFN response in tumors can be attributed to the innate sensing of cytosolic RNA derived from transposable elements (TEs) 46. Considering their role as internal adjuvants eliciting T cell responses against tumors, we hypothesized that differences in TE expression levels could directly influence the immune phenotypic transition from T cell-rich to T cell-excluded. To examine this hypothesis, we conducted a spatial analysis of the expression of 40 TEs that specifically over-expressed in CRC tumor cells (Fig. 4H and Methods) 47. Our analysis identified two CRC-specific TEs—the internal sequence of human endogenous retrovirus subfamily H (HERV-H) and its associated long terminal repeat LTR7Y 48, 49—were significantly downregulated in clone3 compared to clone1 and clone2 (Fig. 4, I and J, fig. S5H, and Methods). The statistical analysis across all slices verified the observations in CRC4 (fig. S5I). This downregulation suggests that tumor cells might suppress HERV-H expression to evade immune recognition in CRC.

To investigate whether the decreased HERV-H expression and attenuated antiviral response observed in clone3 were associated with chr20q-gain, we analyzed TCGA-COAD tumor samples. Notably, tumors with chr20q-gain showed significantly reduced HERV-H expression and diminished viral transcription and type I IFN signaling compared to chr20q-WT (Fig. 4K, fig. S5, J to K, and Methods). We also confirmed the attenuated type I IFN response in tumor cells with chr20q-gain using scRNA-seq data from our cohort (fig. S5, L to O). These results collectively suggest a compelling correlation between the chr20q-gain in tumor cells and the downregulation of HERV-H, which may be implicated in the immune evasion strategies of CRC.

Given our discovery of a strong correlation between chr20q-gain and the consequent downregulation of HERV-H, which results in a T cell-excluded immune phenotype, coupled with the recognized influence of T cells on the clinical response to immune checkpoint blockade (ICB) treatment 50, we hypothesized that chr20q-gain might influence the response to ICB treatment in CRC patients. To examine this, we evaluated survival post-ICB treatment by categorizing patients from Memorial Sloan Kettering Cancer Center (MSKCC) IMPACT (MSK-IMPACT) cohorts (comprising 71 patients) 51, 52 into two groups based on their chr20q status. Our evaluation revealed that patients with chr20q-gain exhibited significantly poorer survival outcomes following ICB treatment (Fig. 4L and fig. S5P). In contrast, this adverse survival trend was not observed in CRC patients who did not receive ICB treatment (TCGA-COAD dataset) (Fig. 4L and fig. S5Q). To rule out the potential impact of tumor mutation burden (TMB), we categorized patients based on both CNV status and TMB. Our analysis indicated that the correlation between chr20q-gain and diminished survival in the CRC ICB cohort remained irrespective of TMB levels (fig. S5R). Taken together, we identified chr20q-gain, a prevalent genomic defect in CRC, as a genomic biomarker for T-cell exclusion, suggesting its utility as an immunogenomic parameter for predicting responsiveness to ICB. Furthermore, our findings lay the groundwork for a novel model in which tumor cells with chr20q-gain orchestrate a T cell-excluded clonalTIME by suppressing HERV-H expression during clonal evolution, thereby potentially evading immune surveillance.

UCASpatial accurately annotates cell subpopulations with spatiotemporal dynamics in the murine wound healing ST dataset

To demonstrate UCASpatial’s utility, we next applied it to another real ST dataset, focusing on murine wound healing–a complex, highly organized multicellular process with distinct spatiotemporal dynamics 53-55. We employed a 2 mm ear punch wound model in MRL/MpJ mice (hereafter denoted as MRL), a strain with super-healing capacity, which fully regenerates the ear punch wound without fibrosis 56-58. As a control, we used non-regenerating strain C57BL/6 mice (hereafter denoted as B6), known to form scar tissue over ear wound edges (Fig. 5A) 56-58. To interrogate the tissue-wide spatiotemporal dynamics of various cell subpopulations during wound healing, we collected ear tissue samples from MRL and B6 mice at the following stages: unwounded (day 0), early inflammatory (day 3), return-to-unwounded (day 7), and regeneration stage (day 15) 54. We performed ST sequencing using the 10x Genomics Visium platform (Fig. 5, B and C). As a representative reference for capturing the cellular diversity present throughout the wound healing process, we utilized scRNA-seq data from day 7 ear tissues from both MRL and B6 mice, classifying them into 23 distinct cell subpopulations (Fig. 5D and fig. S6, A and B).

UCASpatial successfully reconstructed the anatomical distributions of structural cell types, including keratinocytes, endothelial cells, and fibroblasts (Fig. 5, C and E, and fig. S6C). Moreover, UCASpatial resolved the spatiotemporal dynamics of wound-associated immune cells. Specifically, an increase in myeloid cells–neutrophils, basophils, and macrophages–was observed at the early inflammatory stage in both MRL and B6 mice (Fig. 5F). This observation is consistent with the established notion that innate immune cells are the first responders to injury 59-61. At the return-to-unwounded stage, neutrophils and basophils exhibited a reduction in frequency, whereas macrophages maintained a high abundance (Fig. 5F), underscoring their roles in both the early injury response and the subsequent tissue repair and regeneration processes 54, 59, 60, 62, 63.

We experimentally validated the spatial distribution and proportion of six major cell types including keratinocytes, fibroblasts, endothelial cells, neutrophils, macrophages, and basophils, using multiplexed immunohistochemistry (mIHC) (fig. S6D). The anatomical distributions of major cell types detected by mIHC were consistent with UCASpatial’s deconvolution analysis (Fig. 5, C and E, and fig. S6, E to G). Notably, the spatiotemporal dynamics of cell subpopulations quantified by UCASpatial were more closely aligned with the mIHC results than those obtained by other methods (Fig. 5F and fig. S6G). This consistency between the mIHC and UCASpatial’s deconvolution results validates the accuracy of UCASpatial.

To compare the performance of UCASpatial with existing methods, we calculated the PCC between the spatial distribution of cell-subpopulation proportions and marker-gene scores (Supplementary Data 2 and Methods), which allowed us to quantify the performance of deconvolution in the absence of ground truth 16. Among all methods compared, UCASpatial exhibited the highest PCC across cell subpopulations (fig. S6, H and I). Specifically, by comparing UCASpatial with “UCASpatial ablation”, which shares the same framework but lacks the core functions of entropy-based weighting and meta-purity filtering process, we noticed that the improved accuracy of UCASpatial is attributed to the incorporation of entropy weighting and meta-purity filtering (fig. S6, H and I). Notably, existing methods lacked precision in identifying immune cell subpopulations during wound healing. For example, RCTD, stereoscope, and Spotiphy failed to accurately capture the spatiotemporal dynamics of neutrophils, which peaked at the early inflammatory stage in both strains (fig. S6G). Additionally, RCTD, cell2location, CARD, and Spotiphy showed limited ability to detect macrophages at the unwounded stage (fig. S6G).

Further, UCASpatial precisely captured the spatial segregation patterns of macrophage subpopulations. Specifically, MHCII+ macrophages, present at all measured time points, were predominantly located in the wound’s proximal and distal regions in MRL mice, while predominantly located in the wound’s bed and distal regions in B6 mice (Fig. 5, G and H, and fig. S6, J to L). Arg1+ Il1b+ inflammatory macrophages, which appeared transiently in the wound bed at the early inflammatory stage (Fig. 5I), preferentially clustered at the migrating epithelial tongue (Fig. 5G and fig. S6M). In contrast to Arg1+ Il1b+ macrophages, Cd163+ macrophages with enriched functions in endocytosis and lysosome (fig. S6N) predominantly resided along the epidermis layer outside the migrating epithelial tongue across all measured time points post-injury (Fig. 5G and fig. S6M). Cd36+ Gpnmb+ Il1b– macrophages, functionally resembling lipid-associated macrophages implicated in fibrosis 64-66 (fig. S6N), are located in the wound proximal and distal regions at the return-to-unwounded stage (Fig. 5, G and H). Such spatial distribution patterns of four macrophage subpopulations at the return-to-unwounded stage were confirmed by mIHC (Fig. 5I).

We also observed a differential distribution pattern of Cd36+ Gpnmb+ Il1b- macrophages in B6 and MRL mice. This subpopulation was initially located in the wound proximal and distal regions of both strains at the return-to-unwounded stage (Fig. 5, G and H, and fig. S6L). However, at the regeneration stage, Cd36+ Gpnmb+ Il1b- macrophages in B6 mice accumulated within the wound bed, while in MRL mice, they remained distal to the wound bed (Fig. 5H and fig. S6O). This unique migration pattern observed in B6 mice might explain their scarring-healing phenotype. Collectively, UCASpatial reconstructs a fine-grained spatiotemporal cellular landscape in murine wound healing, revealing the unique spatial segregation and temporal dynamics of macrophage subpopulations that serve in different stages of wound healing.

Spatial characterization of molecularly distinct chondrocyte subpopulations linked to regenerative versus scarring-healing phenotypes

After demonstrating our capability to spatially map various cell subpopulations in wounded tissues, we proceeded to investigate the cellular events orchestrating regenerative versus scarring-healing phenotypes in MRL and B6 mice between 7 to 15 days post-injury. We identified a scarring-associated myofibroblast subpopulation Lrrc15+ fibroblasts located within the wound bed at the return-to-unwounded stage in both MRL and B6 mice 67-69 (Fig. S7, A and B). While a reduction in the frequency of Lrrc15+ fibroblasts was observed in both strains at the regeneration stage, these cells were significantly more enriched in B6 mice compared to MRL (Fig. S7, A and B). This differential enrichment may account for the scarring-healing phenotype observed in B6 mice.

We also delineated the spatiotemporal dynamics of keratinocytes and chondrocytes, two major cell populations contributing to regenerative healing 70, 71. We observed a comparable distribution of keratinocytes in MRL and B6 mice within regions of the migrating epithelial tongue at the return-to-unwounded stage (Fig. S7, C and D, and fig. S7E), indicating complete re-epithelialization occurred in both strains. However, at the regeneration stage, while the proportion of keratinocytes remained high within the wound bed of MRL mice, a substantial decrease was observed in B6 mice, supporting the notion that keratinocytes play a crucial role in the regeneration process (Fig. S7, C and D) 72-74.

Spatial mapping of four chondrocyte subpopulations (fig. S7, F and G) during post-injury chondrogenesis revealed distinct regional and temporal patterns. Injury-induced, de-differentiated Spp1+ chondrocytes 75 were acutely present at the wounded bed, with comparable levels in MRL and B6 mice at the early inflammatory stage (Fig. 6A). Prg4+ chondrogenic progenitors, which give rise to hypertrophic chondrocytes 76-78, were transiently observed within the wound bed of B6 mice but persisted in MRL mice from the return-to-unwounded stage till the regeneration stage (Fig. 6, A and B). Consistently, a higher proportion of mature hypertrophic chondrocytes was observed in MRL mice at the regeneration stage (Fig. 6A), suggesting that this progenitor is linked to regenerative healing in MRL mice.

In addition to the Prg4+ chondrogenic progenitor, we noticed an uncharacterized Igfbp5+ chondrocyte subpopulation that expressed Cd34 and Ogn (fig. S7G). Igfbp5+ chondrocyte was closely related to Prg4+ chondrogenic progenitor and was maintained in a less-differentiated state, as evidenced by lower expression levels of Col2a1 and Acan (fig. S7, H to K) 79, 80. Functional analysis of Igfbp5+ chondrocytes showed enriched pathways for extracellular matrix (ECM) organization and ECM-receptor interaction (Fig. S7L), implying potential interaction between Igfbp5+ chondrocytes and ECM-producing cells like fibroblasts and macrophages 81. Intriguingly, the spatiotemporal dynamics of Igfbp5+ chondrocytes differed in MRL and B6 mice. While these cells initially occupied regions proximal to the wound edge in both strains at the return-to-unwounded stage, they continue spreading into the wound bed, replacing the Prg4+ chondrogenic progenitor at the regeneration stage in B6 mice (Fig. 6, A and B). In contrast, Igfbp5+ chondrocytes were spatially restricted in MRL mice, remaining distal to the wound bed at the regeneration stage (Fig. 6, B and C). Therefore, the replacement of Prg4+ chondrogenic progenitors with Igfbp5+ chondrocytes is associated with the scarring healing phenotype in B6 mice.

Igfbp5+ chondrocyte, Cd36+ Gpnmb+ Il1b- macrophage and Fmod+ fibroblast formed a three-cell-type cluster underlying scarring wound healing

Next, we asked whether Igfbp5+ chondrocytes positioned within the wound bed could locally crosstalk with surrounding cells and form a niche that favors scarring healing in B6 mice. To test this, we calculated the PCC between 23 cell subpopulations and identified two chondrocyte-related communities by hierarchical clustering (Methods). The community1, comprising Spp1+ chondrocytes, neutrophils, and Arg1+ Il1b+ macrophages, formed a transient “inflammatory community” within the wound bed at the early inflammatory stage in both MRL and B6 mice (Fig. 6D and fig. S7M). The community2, which we termed the “pro-fibrotic community”, featured a stable three-cell-type cluster involving Igfbp5+ chondrocytes, ECM organization-related Fmod+ fibroblasts (fig. S7N), and Cd36+ Gpnmb+ Il1b- macrophages (Fig. 6, D to F, and fig. S7J). This cluster was particularly prominent in the wound bed at the regeneration stage in B6 mice, which has been validated by mIHC (Fig. 6, G and H, and fig. S7O).

Analysis of functional-enriched pathways across three cell subpopulations within community2 suggested their overlapping functions were responsible for fibrosis, including extracellular matrix organization, collagen fibril organization, and supramolecular fiber organization (fig. S7P and Methods). To determine whether the scarring-healing phenotype observed in B6 mice could be attributed to the formation of this pro-fibrotic community, we conducted cellular-level intervention of the pro-fibrotic community by blocking Fmod+ fibroblast differentiation. Trajectory analysis using the diffusion map algorithm revealed that Fmod+ fibroblasts might originate from Cd34+ fibroblasts (fig. S7Q). During this transition, PPARγ signaling-pathway-related genes were down-regulated, which are known to be associated with fibrosis and regeneration failure (fig. S7R) 82-84. As expected, treating B6 mice with the PPARγ agonist rosiglitazone inhibited the differentiation of Cd34+ fibroblasts into Fmod+ fibroblasts (fig. S7, S and T, Methods). We observed that the three-cell-type cluster in community2 was prominent in the wound bed of control B6 mice at the regeneration stage, whereas this cluster was absent in rosiglitazone-treated B6 mice (fig. S7, U and V). Additionally, rosiglitazone treatment led to a significant reduction in the remaining ear hole area, indicating enhanced wound healing (fig. S7W). Together, these findings underscore the spatial coordination among Igfbp5+ chondrocyte, Cd36+ Gpnmb+ Il1b- macrophage, and Fmod+ fibroblast, forming a three-cell-type cluster that contributes to the scarring-healing phenotype observed in B6 mice.

Disruption of the IL11-IL11RA axis within the pro-fibrotic community enhances wound healing in B6 mice

Next, we examined the potential cross-talks within the pro-fibrotic community. Analysis of specific ligand-receptor interactions revealed a stable interaction between Fmod+ fibroblasts and Cd36+ Gpnmb+ Il1b- macrophages (Fig. 6, I and J), supporting their known partnership 85. Additionally, we observed reciprocal interactions between Igfbp5+ chondrocytes and both Fmod+ fibroblasts and Cd36+ Gpnmb+ Il1b- macrophages. Specifically, Igfbp5+ chondrocytes can crosstalk with Fmod+ fibroblasts via the ligand-receptor pair Timp1-Cd63 (Fig. 6I), the ligand-receptor pair for the induction of collagen synthesis during fibrosis 86, 87. Further, Igfbp5+ chondrocytes can potentially interplay with Cd36+ Gpnmb+ Il1b- macrophages through the ligand-receptor pair of C4-C3ar1 and Csf1-Csf1r (Fig. 6I), signals required for macrophage survival and polarization 88, 89. Cd36+ Gpnmb+ Il1b- macrophages and Fmod+ fibroblasts exhibited interaction with Igfbp5+ chondrocytes via Il11-Il11ra1, the ligand-receptor pair that could promote fibrosis (Fig. 6I) 90. Besides, Cd36+ Gpnmb+ Il1b- macrophages might also promote the proliferation of Igfbp5+ chondrocytes and Fmod+ fibroblasts via the interaction of Grn-Tnfrsf1a (Fig. 6I) 91.

Among these ligand-receptor interactions, we identified an unexplored axis involving the Il11-Il11ra1 signaling pathway in the pro-fibrotic community (Fig. 6I). IL-11 is a member of the IL-6 cytokine family and plays a pivotal role in driving organ fibrosis and dysfunction 92-94. Given that IL11RA is expressed in a wide array of cell types, including stromal cells, pericytes, and epithelial cells, IL-11 signaling can trigger diverse cellular responses, leading to varied outcomes depending on the affected cell type 92. It remains unknown which specific cell types receive IL-11 signals during wound healing and whether this subsequently affects wound healing.

We noticed that Il11ra1 expression was dramatically up-regulated in Cd36+ Gpnmb+ Il1b- macrophages, Igfbp5+ chondrocytes, and Fmod+ fibroblasts of B6 mice compared to MRL mice (Fig. 6K and fig. S7X). This strain-specific difference suggests that the IL11-IL11RA axis might be responsible for the formation of the pro-fibrotic community in B6 mice. To test this, we blocked IL11-IL11RA signaling using an anti-IL11RA neutralizing antibody (anti-IL11RA Ab) in B6 mice. We found that the formation of the pro-fibrotic community was prevented in the wound bed in mice receiving anti-IL11RA Ab compared to the control group, indicating that IL11-IL11RA-mediated interactions contribute to the formation of the pro-fibrotic community in B6 mice (Fig. 6L and fig. S7V). Furthermore, neutralization of IL11RA significantly reduced the remaining ear hole area, suggesting that IL-11 signaling limits regeneration in B6 mice (Fig. 6M). Together, these results reveal that the strain-specific difference in wound healing between B6 and MRL mice is driven by the IL11-IL11RA1 axis, underscoring the power of UCASpatial in dissecting the molecular mechanisms underlying complex phenotypes.

Methods
UCASpatial model

In brief, UCASpatial is an entropy-weighted model that annotates the cell-identity (equal to cell-type) compositions on each spatial location by deconvoluting the spatial expression count matrix into a predefined set of reference cell types. The model requires both spatial transcriptomics data and paired or unpaired scRNA-seq data from a given tissue as input. By default, the model takes an untransformed spatial-expression count matrix 𝐏 ∈ ℝ𝑆×𝐺 of genes 𝑔 = {1,..,𝐺} at spots 𝑠 = {1,..,𝑆} as input, as, for example, obtained from the 10x Space Ranger software (10x Visium data). The scRNA-seq data 𝐁 ∈ ℝ𝐺×𝑁 serve as a reference and consist of 𝑁 cells with the same set of 𝐺 genes with 𝐏. All the 𝑁 cells will be classified into 𝐾 cell types which need to be provided by the users.

We denote 𝐗 ∈ ℝ𝑆×𝐾 as the cell-type composition matrix, where each row 𝑥 represents the proportions of the K cell types on spot s. By default, this matrix is estimated using an entropy-weighted non-negative least square model, which allows for accurate and robust annotation in complex situations. (see below).

Meta-purity filtering

The existence of intermediate cells blurs cell-type borders due to their mixed identities from multiple cell types, which will disturb the recognition of cell-identity-specific gene expression profiles (cGEPs). To remove the effect brought by intermediate cells in the reference scRNA-seq data, we designed the meta-purity filtering strategy as an optional preprocessing.

Based on the input scRNA-seq reference matrix 𝐁, we construct meta cells using the function ‘FindClusters’ from the R package ‘Seurat’ with a default resolution of 100. We recommend the resolution parameter by an expectation of an average of 100 cells per meta-cell, which serves as a parameter selection criterion. Using pre-existing cell-type labels, we define the purity of meta-cells by 𝑝𝑢(𝑀𝐶 ) = max 𝑛𝑐𝑡𝑗, where 𝑀𝐶 denotes the ith meta cell, 𝑝𝑢(𝑀𝐶 ) denotes the purity of 𝑀𝐶 , 𝑛 denotes the cell number of cell type j in 𝑀𝐶 , and 𝑛 denotes the cell number of 𝑀𝐶 . The cells contained in meta cells with purity below the given threshold (default 0.95) will be filtered out.

Notably, the meta-purity filtering process does not result in the loss of any predefined cell types (fig. S1, S to U, fig. S3A, and fig. S5B). The meta-purity filtering process has specific designs to prevent the loss of specified cell subsets. We have implemented a stringent criterion for the meta-purity filtering process. If any loss of the specified cell subsets is detected, the process is immediately halted, and displays an error message prompting a review of input parameters. This ensures that there is no unintended loss of any cell types specified by the user.

Downsampling

Another optional preprocessing for referenced scRNA-seq is down sampling. By default, we will randomly downsample the scRNA-seq reference matrix 𝐁 to 100 cells for each cell type based on a given random seed to improve the speed for downstream analysis. Thus, the processed scRNA-seq expression matrix 𝐁̃ ∈ ℝ𝐺×𝑁̃ has been established, in which 𝑁̃ = 100×𝐾 by default.

Optimized non-negative matrix factorization (NMF)

In order to efficiently extract the cGEP matrix 𝐄 ∈ ℝ𝐺×𝐾 , we utilize non-negative matrix factorization (NMF) on the processed scRNA-seq expression matrix 𝐁̃. Based on the provided cell-type labels, we initialize 𝐇 ∈ ℕ𝐾×𝑁̃ as a one-hot matrix representing the identity of each cell, thus guiding it towards a biologically relevant result. This initialization also aims at reducing variability and improving the consistency between runs. We iterated through NMF to estimate the cGEP matrix 𝐄:

𝐁̃ ~ 𝐄×𝐇 (1)

where each element in the matrix 𝐄 and 𝐇 is constrained to be non-negative.

Entropy-based cGEPs weighting

Considering that the indicative capacity for cell types of each gene is unequal, quantifying the indicative capacity of each gene and taking it into account in the decomposition step would benefit from distinguishing the heterogenous cell states with subtle transcriptome differences. Shannon entropy for the expression percentage of each gene across different cell types can be a great quantification for the indicative capacity of each gene. Thus, we formulate Shannon entropy-based weights for each gene 𝑔 of the extracted cGEP matrix 𝐄 to highlight genes indicative of cell identities.

UCASpatial calculates the expression percentage matrix 𝐏𝐂𝐓 ∈ ℝ𝐾×𝐺 based on the processed scRNA-seq expression matrix 𝐁̃:

𝐏𝐂𝐓 = 𝑛𝑢𝑚(𝐁̃ > 0) / 𝑛𝑢𝑚(𝐁̃) (2)

where 𝐏𝐂𝐓 denotes the element of 𝐏𝐂𝐓 representing cell type k for gene g, 𝑛𝑢𝑚(𝐁̃ ) represents the number of cells in cell type k in 𝐁̃, 𝑛𝑢𝑚(𝐁̃ > 0) represents the number of cells in cell type k with gene g expression over 0. The expression percentage matrix 𝐏𝐂𝐓 will be further normalized into 𝐏̃𝐂𝐓 via each element 𝐏𝐂𝐓 being divided by the sum of the corresponding row:

𝐏̃𝐂𝐓 = 𝐏𝐂𝐓 / ∑𝐾 𝐏𝐂𝐓 (3)

Based on the 𝐏̃𝐂𝐓 matrix, UCASpatial further calculates the Shannon entropy for each gene:

Entropy(𝑔) = −∑𝐾 𝐏̃𝐂𝐓 ×log𝐏̃𝐂𝐓 (4)

For better understanding, UCASpatial provides an adjusted entropy using the maximum entropy minus the entropy for each gene g :

adj.ent(𝑔) = 𝐸𝑛𝑡.𝑚𝑎𝑥𝑚−Entropy(𝑔) (5)

where 𝐸𝑛𝑡.𝑚𝑎𝑥𝑚 = sup(Entropy(𝑔)) = −∑𝐾 1 ×log 1 represents the supremum of Entropy(𝑔). The maximum entropy 𝐸𝑛𝑡.𝑚𝑎𝑥𝑚 is only related to the number of cell types K. The adjusted entropy could excellently represent the expression specificity for genes, with a higher adjusted entropy representing a higher indicative capacity for cell types of a gene.

Considering that the adjusted entropy can only evaluate the expression specificity of each gene, we thus took not only the adjusted entropy but also the expression percentage and expression level into account, designing the entropy-based weight of each gene as:

𝜔(𝑔) = max ( 2×adj.ent(𝑔) / 𝐸𝑛𝑡.𝑚𝑎𝑥𝑚 ×2×𝐏̃𝐂𝐓 ×avg_logFC(𝑖,𝑔)) (6)

𝜔(𝑔) = { 𝜔(𝑔), 𝜔(𝑔) ≤ 1; 1, 𝜔(𝑔) > 1 } (7)

where avg_logFC(𝑖,𝑔) denotes the average of log fold-change expression for gene g in cell type i, 𝜔(𝑔) denotes the entropy-based weight of gene g which is constrained to be non-negative.

UCASpatial also combined common marker finding methods including the FindAllMarkers function in Seurat (version 3.1.4) and cosine similarity method 96 (version 0.9.0), for optional gene filtering in cGEPs, which can further improve the indicative capacity of genes in cGEPs.

Entropy weighted non-negative least square model

Taking the untransformed spatial expression count matrix 𝐏 as input, this step aims to decompose the RNA counts of each spot into referenced cell type composition 𝐗 ∈ ℝ𝑆×𝐾 based on 𝐄. Unit variance normalization by gene is performed in 𝐏 in order to standardize discretized gene expression levels 5, 97. Commonly, this step will utilize the non-negative least square (NNLS) model for each spot s:

𝐏 = 𝐗 ×𝐄+𝜖 (8)

where 𝐏 denotes the spatial expression count vector of spot s, 𝜖 ∈ ℝ𝑆×1 represents the residual error of spot s which is restricted to being non-negative. By minimizing the loss function for each spot s:

min𝐿(𝜖) = min‖𝜖 ‖2 = min∑𝐾 ∑𝐺 (𝐏 −𝐗 𝐄 )^2 (9)

where 𝐗 denotes the cell type composition of spot s, with each element restricted to be non-negative; 𝐏 denotes the spatial expression count of gene j in spot s; 𝜖 denotes the residual error between the deconvoluted expression vector 𝐗 𝐄 and true spatial expression 𝐏 . By minimizing the loss function 𝐿(𝜖) = ‖𝜖 ‖2, we can obtain the 𝐗 which can best mimic the real cell type composition of spot s under the given cell identities.

To further take the entropy-based weight into account, we extended the common NNLS into an entropy-Weighted Non-Negative Least Squares (WNNLS) model. We first extend the weight vector 𝜔(𝑔) to a diagonal matrix 𝛀 ∈ ℝ𝐺×𝐺:

𝛀 = diag(𝜔(1), ..., 𝜔(𝐺)) (10)

the corresponding WNNLS model has been extended into:

𝐏 √𝛀 = 𝐗 𝐄√𝛀 (11)

The loss function of the WNNLS model weighted by the diagonal matrix 𝛀 is:

Min𝐿̃(𝜖) = min‖𝜖̃ ‖2 = min‖𝐏 √𝛀−𝐗 𝐄√𝛀‖^2 = min∑𝐾 ∑𝐺 (𝐏 √𝛀 −𝐗 𝐄 √𝛀 )^2 (12)

where each element in 𝐗 is restrained to be non-negative. Each row of the cell-type composition spectrum 𝐗 is normalized into [0,1] by 𝐗 = 𝐗 / ∑𝐾 𝐗 .

By minimizing the weighted loss function min𝐿̃(𝜖) = min‖𝜖̃ ‖2 , we can highlight the contribution of genes indicative of cell identities, and finally calculate the cell-type composition spectrum 𝐗 in a more biological way.

Referring to SPOTlight, by selecting all cells in 𝐁̃ from the same cell type k, and computing the median of each column of 𝐄 for a consensus cell type-specific topic signature, we use 𝐇 to further correct 𝐄.

Additionally, while UCASpatial is implemented in R, its deconvolution algorithm NMF is implemented with an efficient C++ code that is linked back to the main functions of UCASpatial through Rcpp, ensuring scalable computation.

Finally, considering that the number of cells in each spot is generally within a certain range, for example, 5-20 cells for 10x Visium, we believe that a proportion less than the threshold will be noise in the calculation process. By default, we use a less stringent filter of 0.01 to remove potential noise in the cell-type composition 𝐗.

Optional other weighting strategies

In addition to an entropy-based weighting strategy, we implemented alternative weighting strategies to accommodate diverse user requirements and analytical contexts—including mutual information (MI) weighting and variance-based weighting—are provided and described as follows.

Variance-based Weighting: For each gene, mean expression levels were computed across predefined cellular subpopulations (clusters) using normalized scRNA-seq data from the selected assay and data slot. The variance of these cluster-wise mean expressions served as a specificity metric, reflecting the extent of differential expression between subpopulations. These variances were linearly normalized to the [0,1] range to obtain variance-based weights. To integrate gene expression magnitude and percentage, these weights were combined multiplicatively with average log fold change and cluster-wise expression frequency metrics. The final combined weight was bounded to a maximum value of 1.

Mutual Information-based Weighting: To capture nonlinear or complex dependencies between gene expression and subpopulation identity, gene expression was binarized (expressed vs. non-expressed) across cells. Mutual information between this binary expression state and cluster labels was computed using empirical entropy-based estimators, quantifying the reduction in uncertainty of cluster identity given gene expression. Mutual information scores were rescaled to [0,1] and combined with log fold change and expression percentage as above. The resulting weight was capped at 1.

Notably, the weighting strategy can be conveniently adjusted via the weight.strategy parameter interface provided in UCASpatial.

Implementation of the reference diagnosis tool

To systematically assess the quality and specificity of the scRNA-seq reference prior to deconvolution, and to address potential biases arising from dissociation or cell type definition, we implemented a reference diagnosis tool within the UCASpatial framework. This tool evaluates whether the user-provided reference cell subpopulations possess distinct, recoverable transcriptional features under the algorithm's feature extraction mechanism.

Specifically, the tool leverages the correspondence between the extracted latent topics (cGEPs) and the user-annotated cell subpopulations. We visualize this relationship using a dot plot matrix, where the x-axis represents the reference cell subpopulations and the y-axis represents the extracted NMF topics. The size and color intensity of each dot reflect the normalized weight or contribution of a specific topic to a cell subpopulation.

This tool provides a practical and systematic guideline for users: if the diagnosis indicates a poor mapping, the reference requires refinement (e.g., re-clustering, removing low-quality cells, or acquiring an alternative reference).

Simulation strategy and evaluation metrics

Simulation strategy

We generated simulated spatial transcriptomic data sets from a reference CRC scRNA-seq data set (fig. S1, A to D, see in the Human specimens (in-house cohort) part). To simulate different situations of tumor microenvironment, we generated three different regions including tumor core (TC), tumor margin (TM), and tumor stroma (TS) 29, 41. Each region has five replicates with each replicate containing 𝑁 = 50 spots, generated based on random seeds (19980205, 19990911, 19960330, 20000502, 20230509). The number of cells in each spot is restricted to less than 20. To mimic the true situation as much as possible, we also imposed restrictions on the types of cells present in different regions, so that spots from different regions have different cell type co-localization patterns. Detailed restrictions for regional simulation are listed in Supplementary Table 1. To better mimic the real capture, if the synthetic spot had over 25,000 UMI counts, it would be randomly down-sampled into 20,000 NUMI counts.

We applied different methods (see in the Compared methods part) to deconvolute the above datasets. In each analysis, we provided identical simulated spatial data and reference scRNA-seq data as input for all methods. To ensure an equitable assessment of algorithm performance, we adopted the identical threshold of 0.05 (established on spots containing up to 20 cells) as UCASpatial for algorithms lacking default filtering thresholds, aiming to eliminate noise in accessing the ability to determine whether individual cell types exist at a spot (i.e., spots with non-zero cell abundance or not).

Construct complexity hierarchy

To evaluate UCASpatial in microenvironments with different complexity, we rank the complexity hierarchy of these simulated spots by calculating their cell-type composition entropy:

Entropy(𝑠) = −∑𝐾 𝐗 ×log𝐗 , where 𝐗 denotes the proportion of cell type i in spot s. We classify them into three groups (low, medium, and high), with 𝑁 = 250 in each group, where higher entropy indicates a higher level of complexity.

Evaluation metrics

In each simulation replicate, we calculated the true cell-subpopulation proportions on each spatial location as the number of cells in each cell subpopulation divided by the total number of cells in the location. 𝑇 and 𝑃 represent the true and predicted proportion of cell subpopulation k, respectively, in the sth spot through the number of total cell subpopulations K. To evaluate the performance of the methods comprehensively, the main benchmark metrics used were PCC, RMSE, and F1 score. The definitions of these metrics are provided below.

PCC. The cell subpopulation-level PCC value was calculated between the predicted proportion of the specific cell subpopulation 𝑃 and the ground truth proportion of the specific cell subpopulation 𝑇 using the following equation:

PCC (𝑃 ,𝑇 ) = (E[𝑃 ,𝑇 ]−E[𝑃 ]E[𝑇 ]) / (√E[𝑃 2]−(E[𝑃 ])2√E[𝑇 2]−(E[𝑇 ])2) (13)

The spot level PCC value was calculated between the predicted spot composition of 𝑃 and the ground truth spot composition 𝑇 using the following equation:

PCC (𝑃,𝑇 ) = (E[𝑃,𝑇 ]−E[𝑃]E[𝑇 ]) / (√E[𝑃2]−(E[𝑃])2√E[𝑇 2]−(E[𝑇 ])2) (14)

RMSE. The cell subpopulation-level RMSE was calculated between the predicted proportion of the specific cell subpopulation 𝑃 and the ground truth proportion of the specific cell subpopulation 𝑇 𝑘 using the following equation:

RMSE (𝑃 ,𝑇 ) = √ (1/𝑆) ∑(𝑇 −𝑃 )^2 (15)

The spot level RMSE value was calculated between the predicted spot composition of 𝑃 and the ground truth spot composition 𝑇 using the following equation:

RMSE (𝑃,𝑇 ) = √ (1/𝐾) ∑(𝑇 −𝑃 )^2 (16)

F1 score. F1 score was selected to evaluate the ability of alternative methods to determine the presence of individual cell types at specific locations, i.e., with non-zero cell abundances, based on the confusion matrix. A confusion matrix is a table used in the field of machine learning and statistics to assess the performance of a classification model. It summarizes the performance of a classification algorithm by showing the counts of true positive (TP), true negative (TN), false positive (FP), and false negative (FN) predictions for each class in a dataset.

F1 score = (2×(𝑇𝑃/(𝑇𝑃+𝐹𝑃))×(𝑇𝑃/(𝑇𝑃+𝐹𝑁))) / ((𝑇𝑃/(𝑇𝑃+𝐹𝑃))+(𝑇𝑃/(𝑇𝑃+𝐹𝑁))) (17)

Compared methods

We compared UCASpatial with five deconvolution methods: (1) SPOTlight (version 0.1.7), (2) RCTD (version 1.2.0) (3) cell2location (version 0.8a0), (4) CARD (version 1.0) and (5) stereoscope (v 0.3.1). (6) Spotiphy (version 0.1.3). Results of the stereoscope were computed with parameters “-sce 25000 -n 5000 -ste 25000 -stb 500 -scb 500 -gp”. Results of the RCTD were use the ‘full model’. For other methods, we followed the tutorial on the corresponding GitHub pages and used the recommended default parameter settings for deconvolution analysis.

Computer platform for evaluating runtime and memory usage

We ran CPU tests of the 7 methods on a computer cluster with an identical CPU platform (2.1 GHz, 60.5 MB L3 cache, 44 CPU cores) and 384 GB memory (DDR4 2,666 MHz). The GPU tests for cell2location and stereoscope were performed on the same computer with an additional GPU (GeForce RTX 4090 with 24 GB memory). To assess the impact of various data attributes (including the number of cells in scRNA-seq data and the number of spots in spatial data) on the computing resources consumed by those 7 methods, we constructed simulated spatial datasets with: 1) varying number of spots (50, 500, 5000, and 50000), 12500 cells, and 20 cell types; 2) varying numbers of cells (approximately 1000, 2500, 5000, 10000, 15000, 25000, and 35000), 20 cell types, and 5000 spots. For the memory evaluation, we used the simulated dataset with 50000 spots ST and about 12500 cells, and 20 cell types. Then, we evaluated the impacts of these data attributes on computing resources consumed.

Analysis of human CRC 10x Visium HD data

To test UCASpatial’s ability on the 10x Visium HD platform, we applied UCASpatial to 10x Visium HD data from a human CRC sample 31. We used the recommended bin size of 8 μm × 8 μm for this analysis. With the scRNA-seq data shown in Fig. 3B as a reference, we spatially annotated 34 cell subpopulations using the hyper-parameter setting ‘meta.filter = F’ (fig. S2A). For comparison, we also performed RCTD using the ‘full’ model. Other aforementioned methods did not finish the task for the 8μm bin annotation, possibly due to high memory demands, long run times, or software version incompatibility with 10x Visium HD. The high-resolution images of the deconvolution results have been uploaded to GitHub (https://github.com/BIGHanLab/UCASpatial) for better visibility.

In our testing environment (2.1 GHz, 60.5 MB L3 cache, 44 CPU cores), UCASpatial completed the task in just 1.97 hours, while RCTD required 2.27 days for complete annotation. Furthermore, we evaluated UCASpatial on 10x Visium HD data with 2 μm × 2 μm bin size, which was processed in 2.14 hours, demonstrating its scalability to high-resolution data.

Analysis of mouse hippocampus Slide-seq V2 data

To test UCASpatial’s ability in Slide-seq V2 platform, we used UCASpatial (meta.resolution = 200, meta.purity = 0.98, cos.filter.threshold = 0.1, ent.filter.threshold = 0.2, logfc.threshold = 0.1, min.diff.pct = 0.05, meta.assay = ‘RNA’) to spatially annotate cell subpopulations in Slide-seq V2 data of the mouse hippocampus using published scRNA-seq datasets of mouse hippocampus10, 98. It is noteworthy that for CARD, we utilized a downsampled scRNA-seq reference dataset consistent with the original CARD paper 12 to faithfully reproduce its results, whereas other methods were run using the full single-cell reference. Similar to the strategy in the CARD paper, we didn’t apply cell2location on this dataset due to its heavy computational burden. Results of the stereoscope were computed with parameters “-sce 25000 -n 5000 -ste 25000 -stb 500 -scb 500 –gp”. Results of the RCTD were used with the ‘full model’. For other methods, we followed the tutorial on the corresponding GitHub pages and used the recommended default parameter settings for deconvolution analysis. The visualization results of each method can be found in fig. S2H.

Simulated data generation and method benchmarking based on 10x Visium HD data

Public 10x Visium HD data from human colorectal cancer (CRC) samples were obtained31. Cell segmentation and spatial alignment were performed using Space Ranger software (version 4.0.1), and reads were aligned to the hg38 reference genome. Cells were annotated according to their transcriptomic profiles (fig. S2, J and K, and Fig. 3A). To approximate the capture area of a standard Visium spot (55 µm in diameter), we overlaid a grid of 49 µm × 49 µm square bins onto the coordinate system of each tissue section. For each square bin, we aggregated the transcriptomic counts (UMI counts) from all segmented cells whose centroids lay within its boundaries, thereby generating a composite gene expression profile representing a single “pseudo-spot” (fig. S2I).

To evaluate the performance of each deconvolution method on these pseudo‑spots, we employed UCASpatial (with meta.filter = FALSE) to spatially annotate cell subpopulations within each pseudo‑spot. We did not apply Cell2location to this dataset due to its substantial computational demands. Results from Stereoscope were obtained using the parameters “-sce 500 -n 5000 -ste 500 -stb 500 -scb 500 –gp.” For RCTD, we used the “full model.” For all other methods, we followed the tutorials available on the respective GitHub repositories and applied the recommended default parameter settings for deconvolution analysis.

Detailed deconvolution results of each method are publicly available in the Zenodo repository (see Data Availability).

Experimental and Analytical Approach to Human CRC

Human specimens (in-house cohort)

This study complies with all relevant ethical regulations and was approved by the Institutional Ethics Committee of Beijing Cancer Hospital (approval number: 2023YJZ16). All patients were enrolled and pathologically diagnosed with colon cancer at Peking University Cancer Hospital & Institute (Supplementary Table 2). Written informed consent was obtained from all patients. Clinical and pathological stages were determined according to the guidelines of the AJCC version 8. Detailed clinical information was summarized in Supplementary Table 2. Fresh tumor tissues and patient-matched normal adjacent colon tissues were collected from surgical resection specimens, with no prior exposure to chemotherapy or radiation therapy. The resected samples were stored in RPMI 1640 (Gibco) supplemented with HEPES (Thermo Fisher, 15630080), Penicillin-Streptomycin (Thermo Fisher, 15140122), Amphotericin B (Sigma, S1636), and Gentamicin Sulfate (Gibco, 15710-064) and transported at 4°C immediately. As for CRC1, CRC2 and CRC3, a small portion of each tissue was frozen and embedded in OCT (Sakura Tissue-Tek), while the remainder of the same tissue was dissociated into a single-cell suspension.

Single-cell RNA-seq for patient samples

Fresh human tumor specimens and patient-matched adjacent normal tissues were washed twice with DPBS (gibco) and then mechanically disrupted into approximately 1 mm3 pieces in the RPMI-1640 (gibco) and further enzymatically digested with Human Tumor Dissociation Kit (Miltenyi Biotec, 130-095-929) in C-tube in combination with gentleMACS™ Dissociator (Miltenyi Biotec, 130-093-235) presetting program h_tumor_01, h_tumor_02 and h_tumor_03 sequentially according to the manufacturer’s instruction. Between each program, samples were incubated for 30 minutes at 37℃ under continuous rotation using a rotator. The dissociated cells were subsequently passed through a 70-μm cell strainer (Biologix, 15-1070) and centrifuged at 300g for 7 min. After the supernatant was removed, the cell pellet was washed with DPBS. The pelleted cells were resuspended in a cell cryopreservation medium (FBS with 10% DMSO) at a concentration of 1×106 cells/mL and further cryopreserved in liquid nitrogen for long-term storage. On the day of scRNA-seq droplet capture, cells were thawed rapidly in a 37℃-water bath and washed off DMSO in DPBS. Palleted cells were resuspended in a sorting buffer (PBS with 2% FBS and 1 mM EDTA). Cells were stained with CD45-FITC (BioLegend Cat# 368508, clone 2D1, RRID: AB_2566368, 100ug/ml, 1:500) for 20 minutes at 4℃ in the dark. 7-Amino-Actinomycin D (7-AAD, Sigma-Aldrich Cat# 129935) was freshly added at a final concentration of 1 μg/mL prior to FACS sorting performed on the BD FACS Aria III instrument. Based on FACS analysis, 7-AAD-CD45+ cells, and 7-AAD-CD45- cells were sorted into 1.5ml tubes (ExCell Bio, CS015-0041) and counted manually under the microscope (Leica), respectively. Libraries were generated using Chromium Next GEM Single Cell 5’ Library & Gel Bead Kit v1.1 (10x Genomics, PN-1000165) according to the manufacturer’s protocol. Purified libraries were sequenced by the Illumina NovaSeq 6000 platform with 150-bp paired-end reads.

10x Visium spatial transcriptomics

CRC samples were frozen and embedded in OCT (Sakura Tissue-Tek). All tissue blocks were sectioned using the Leica CM1950 cryostat and were cut at 10 μm. Samples were selected on the basis of morphology, orientation (based on H&E staining), and RNA integrity number (RIN) that was obtained using Agilent 2100 Bioanalyzer. The RIN of selected samples ranged from 4.5 to 8.3. The sections were treated as recommended by 10x Genomics and the optimization procedure showed the optimal permeabilization time determined as 12 min for CRC samples using the Visium Spatial Tissue Optimization Kit (10x Genomics, PN-1000193). Spatial gene expression slides were equilibrated inside the cryostat at -20℃ before sectioning and all of the reagents used for spatial transcriptomics were from the Visium Spatial Gene Expression Kit (10x Genomics, PN-1000187). According to the manufacturer’s instruction, the tissue sections were fixed in chilled methanol at -20℃ for 30 min and stained with H&E. Brightfield histological images were performed on the Olympus Microscopes BX53 at ×20 magnification and then taken by cellSens Dimension (version 1.18) software. Images were exported as tiled tiffs for analysis. After imaging, tissue sections were permeabilized at 37℃ for optimal permeabilization time of release of mRNA as mentioned above. Following permeabilization, mRNA was reverse-transcribed and underwent second-strand synthesis, denaturation, cDNA amplification, and cleanup by SPRIselect beads sequentially (Visium, 10x Genomics, User Guide CG000239). cDNA amplification cycles were determined by quantitative PCR (qPCR). Library quality was evaluated using TapeStation (Agilent Technologies). Sequencing was performed on the Illumina HiSeq X Ten platform.

Single-cell RNA-seq data processing and annotation

In-house single-cell RNA-seq data from 10x Genomics were processed using Cell Ranger (version 6.0.1), aligning to the hg38 reference genome. Low-quality cells were filtered out if they expressed fewer than 500 or more than 7,000 genes. Cells with mitochondrial gene expression exceeding 20% were also excluded. Batch effects of in-house and public scRNA-seq datasets were corrected using the integration workflow in the Seurat R package (version 3.1.4) 99.

To identify and refine cell subpopulations, we employed a multi-step approach. Initially, we employed log-normalization with a scaling factor of 10,000 to adjust for variations in library size among cells. We then identified top 2000 highly variable genes utilizing the Seurat package and conducted Principal component analysis (PCA) on these genes, retaining the top 20 principal components that determined by the elbow method on a Scree plot. Next, we proceeded with an iterative clustering approach to delineate clusters within the single-cell data. Generally, this iterative process initially identified clusters using the Leiden algorithm, followed by a merging of clusters from cell types that appeared to be identical. In the context of certain cell types characterized by heightened heterogeneity, such as T&NK, myeloid, and epithelial cells, we performed cell-type-specific re-integration based on top 2000 variable features to capture more signature genes that can represent the heterogeneity of subpopulations. Then we performed sub-clustering to further identify subpopulations using Leiden clustering. This iterative clustering process allowed us to refine the initial clustering, enabling the splitting of initial clusters that contained multiple subpopulations into more specific lower-level subpopulations. One T cell and myeloid cluster with few specific markers but exhibiting high expression of cell cycle-related genes was removed, as these genes could not be fully regressed out after cell cycle score scaling. In total, we defined 34 cell types in CRC, and this CRC scRNA-seq reference was used throughout the manuscript to maintain consistency in presenting spatial mapping results (Fig. 3B and fig. S3, B to D).

To infer tumor clonal diversity from scRNA-seq data, CNVs were evaluated using the InferCNV v1.2.2 R package (https://github.com/broadinstitute/inferCNV). CNVs of epithelial cells were estimated using non-epithelial cells as a reference. We used a Scutoff of 0.1 for the minimum average read counts per gene among reference cells, clustered each group of cells separately, and denoised our output (fig. S5L).

Spatial transcriptomics data processing

Raw sequencing reads of spatial transcriptomics were quality-checked and mapped by Space Ranger v1.3.0. The gene-spot matrices generated after ST data processing from ST and Visium samples were analyzed with the Seurat package in R. The processed data were used for downstream analysis of UCASpatial (meta.resolution = 300, meta.purity = 0.9, other parameters using default) and other methods (using default parameters). Using the threshold of 3% to filter the noise after decomposition by each method.

Gene signature score calculated in ST

Selected the top 20 genes of each cell subpopulation provided in Supplementary Data 1 or Supplementary Data 2 as the signature gene list for the murine wound healing ST dataset and human CRC ST dataset, respectively. Used the “AddModuleScore” method in the Seurat R package to calculate the signature score of each cell subpopulation in the ST dataset using the signature gene list.

Cellular distance evaluation in spatial transcriptomics

In order to quantify the distance relationships between different cell subpopulations in spatial transcriptomics data, we developed an algorithm based on the deconvolution results of UCASpatial. To calculate the spatial distance between cell subpopulation A and cell subpopulation B in spatial transcriptomics data, we first conduct a breadth-first search for each spot containing cell subpopulation A to identify the nearest spot containing cell subpopulation B, and record the distance as well as the proportions of subpopulation A and B within their respective spots. Then we calculate a weighted average on each shortest distance, taking into account the proportions of subpopulation A and subpopulation B within the spots simultaneously. This yields the distribution of spatial distances between the two subpopulations.

Multi-color immunohistochemistry (mIHC)

mIHC was performed manually using the Opal 7-Color Manual IHC Kit (AKOYA BIOSCIENCES, Cat# NEL81 1001KT), following the manufacturer’s protocol. Fresh human tumor specimens were excised and fixed in 10% formalin for 24 hours at 4°C. The fixed tissues were then embedded in paraffin and sectioned at 5 µm thickness for histological evaluation. Tissue sections were baked in the oven at 65°C for 1 hour, dewaxed with xylene, and rehydrated through a graded series of alcohol concentrations (100%, 95%, 70%).

Endogenous Peroxidase Blocking Buffer (Beyotime Cat# P0100A) is used in the method to quench endogenous peroxidase activity. Microwave treatment is used for antigen retrieval and removing antibodies after a target has been detected. The slide was boiled in the microwave oven for 45 seconds at 100% power and then 15 min at 20% power, followed by blocking for 10 min at room temperature. Primary antibodies against protein CD14 (Proteintech Cat# 60253-1-Ig, clone 2C1D9, RRID:AB_2881374; 1:200), C1qC (Proteintech Cat# 16889-1-AP, polyclonal, RRID:AB_2878328; 1:100), IGF1 (Proteintech Cat# 28530-1-AP, polyclonal, RRID:AB_2881164; 1:100), SPP1 (Proteintech Cat# 22952-1-AP, polyclonal, RRID:AB_2783651; 1:200), PDPN (Proteintech Cat# 11629-1-AP, polyclonal, RRID:AB_2162067; 1:100), ANXA1 (Abcam Cat# ab214486, clone EPR19342, RRID:AB_2890907; 1:4000), CXCL13 (Abcam Cat# ab246518, clone EPR23400-92, RRID:AB_2894862; 1:1000), PD-1 (Abcam Cat# ab137132, clone EPR4877(2), RRID:AB_2894867; 1:500), CD4 (Abcam Cat# ab133616, clone EPR6855, RRID:AB_2750883; 1:100), FOXP3 (Abcam Cat# ab20034, clone 236A/E7, RRID:AB_445284, 1:1000), IL-17a (Proteintech Cat# 26163-1-AP, polyclonal, RRID:AB_2880409, 1:1000), HLA-DR (Abcam Cat# ab92511, clone EPR3692, RRID:AB_10563656, 1:400), CD20 (Abcam Cat# ab78237, clone EP459Y, RRID:AB_1640323, 1:250), CLEC9A (Proteintech Cat# 55451-1-AP, polyclonal, RRID:AB_2881338; 1:200), CD1c (Abcam Cat# ab246520, clone EPR23189-196, RRID:AB_2894854; 1:1000), CD208 (Abcam Cat# ab271053, clone EPR24265-8, RRID:AB_2943181; 1:500), GZMKE (Abcam Cat# ab282703, clone EPR24601-164, RRID:AB_3107073; 1:500), CD8 (Abcam Cat# ab17147, clone C8/144B, RRID:AB_443686; 1:100), CD103 (Abcam Cat# ab224202, clone EPR22590-27, RRID:AB_2936238; 1:800), CTHRC1 (Abcam Cat# ab256458, clone EPR22851-145, RRID:AB_3096210; 1:400), and TCF1 (Cell Signaling Technology Cat# 2203T, clone C63D9; 1:200) were sequentially applied at 4°C overnight or at room temperature for 1 hour.

The secondary antibody Polymer HRP Ms + Rb (Cat# ARH1001EA) was applied and incubated at room temperature for 10 minutes. Signal amplification was performed using the Opal Fluorophore Working Solution on the slides at room temperature for 10 minutes. Finally, nuclei were stained with DAPI (Invitrogen Cat# D1306, RRID: AB_2629482, 5 mg/mL, 1:2500). Multi-color IHC imaging was collected using the Vectra Polaris Imaging System (PerkinElmer) 100 at 20x magnification and integrated using the InForm Advanced Image Analysis Software (PerkinElmer, version 2.4.1).

Cell segmentation and quantification was carried out using QuPath 101. In short, DAPI stains were used for cell detection with the following parameters: requestedPixelSize 1 μm, backgroundRadius 10 μm, sigma 0.5 μm, minArea 10 μm, maxArea 400 μm, threshold 0.1, cellExpansion 1.5μm. Cell subpopulations were classified based on the combinatorial expression of key protein markers (fig. S4C). The distances between CD8+ T cell subpopulations and the tumor region were calculated using the “Distance to annotation 2D” function in QuPath. For the pseudo-spot visualization, we first divide the mIHC data into densely distributed bins, each with a width and length of approximately 50 μm. Then, we calculate the proportion of each cell subpopulation in the pseudo-spots by dividing the number of specified cell subpopulations by the total cell count defined by DAPI.

Tumor clonal inference and clonalTIME quantification using UCASpatial-clonalTIME

To infer tumor clonal diversity within each spatial transcriptomics slice, CNVs were evaluated by ordering genes according to their chromosomal location and applying a moving average to the relative expression values. Genes were arranged according to their genomic localization using the SpatialInferCNV R package (version 0.1.0) 33. CNVs in high epithelial spots (overall epithelial proportion > 40%) were estimated using reference spots containing few or no malignant cells selected from UCASpatial annotated remaining spots (overall epithelial proportion < 40%). Considering the non-dispersed spatial distribution of epithelial cells, scattered high epithelial spots (neighboring same-type spots < 1, step = 1) were considered noise and defined as remaining regions. The spot-resolution CNV profiles for each spatial transcriptomics slice were then estimated using the SpatialInferCNV workflow (cutoff = 0.1, cluster_by_groups = FALSE, denoise = TRUE, HMM = FALSE), excluding chrM, chrX and chrY genes. To further define tumor regions and infer tumor clonal diversity based on CNV data, we constructed dendrogram trees to determine clonal groupings of each spot. Dendrogram nodes with at least one high-probability CNV region were defined as tumor clones and assigned a specific identity. Nodes without high-probability CNVs were assigned to normal epithelial or remaining regions based on their histological structure. The average cell subpopulation composition across all spots within a clonal region was defined as the cell subpopulation composition of that clonalTIME (Fig. 3H). To perform a comprehensive comparison of individual CNVs across clones, we initially identified the CNV locations and statuses of each clone using an HMM model in SpatialInferCNV (Fig. 4, A to C, and fig. S5A). CNV regions located on the same chromosome with similar variation patterns in the TCGA-COAD dataset were merged if their PCC (log2FC_CNV) > 0.8, where log2FC_CNV denotes the log2 fold change of the copy number of the tumor over the copy number of normal. These CNVs were then separated based on the location of the centromere (Supplementary Table 4).

ClonalTIME phenotyping

To classify the clonalTIME phenotype of each tumor clone, the scaled proportion of the T cell subpopulation of each clone was clustered using the K-means function (centers = 2). ClonalTIMEs were classified as T cell-riched or T cell-excluded based on the proportion of T cell subpopulation within each clustering group (Fig. 3I and fig. S4, J and K). To determine whether a given cell subpopulation was enriched in a specific clonalTIME phenotype, we calculated the ratio of observed cell number to expected cell number (Ro/e) for T cell-rich clonalTIME, where the expected proportion of cell subpopulation was determined from the chi-squared test. We assumed that a subpopulation was enriched in T cell-rich clonalTIME if Ro/e > 1.2 and enriched in T cell-excluded clonalTIME if Ro/e < 0.8.

Multi-cellular community analysis

To identify cell subpopulation co-localization patterns, we first adjusted the cellular composition of each spot by averaging the proportions of cell subpopulations within that spot and its immediate neighbors (step = 1). Subsequently, we employed the Pearson method to calculate the correlations among the adjusted cell subpopulations across all spots within the tumor region. To identify distinct multi-cellular communities, we applied the ward.D2 hierarchical clustering algorithm to the correlation matrix. This analysis yielded three multicellular communities in the human CRC part, and two multicellular communities in the murine wound healing part (Fig. 3J and Fig. 6D). The functionality of each multicellular community was characterized based on the shared enriched pathways exhibited by its constituent cell subpopulations. Shared enriched pathways were defined as those that were statistically significantly enriched in at least half of the subpopulations within each community, with a minimum of one subpopulation per community contributing to this enrichment. To address the issue of redundancy in pathway databases, the identified shared enriched pathways were subjected to hierarchical clustering using the pairwise Pearson correlation distance calculated across the gene space. Here, each pathway was represented by its associated differentially expressed genes. Pathway clusters were named manually.

Multi-cellular community score

The proportions of each cell subpopulation were first normalized by the proportion of non-epithelial cells, ensuring the sum of non-epithelial cell proportions of a spot equals one. To enable inter-cell type comparability, we scaled the proportion of each cell subpopulation using the 99.5th percentile across all tumor spots, resulting in a score for each cell subpopulation. The average score of all cell subpopulations within a community was then calculated as the community score.

Spatial quantification of TE expression

We quantify the expression of transposable elements (TEs) using spatial transcriptomic data based on a previously reported algorithm 47. This algorithm allocates TE reads to TE metagenes based on the TE type-specific sequence (Fig. 4H). In brief, the genes/TEs and cell count matrix were generated according to default parameters and normalized using the NormalizeData function. Clonal-specific TE expression programs were identified by comparing each tumor clone region using the FindAllMarkers function. The list of CRC-specifically upregulated TEs was downloaded from previous literature 48. The HERV-H expression level was calculated as the average expression of previously reported HERV-H elements (HERVH-int, LTR7, LTR7B, and LTR7Y) (Fig. 4J).

Differential expression analysis using spatial transcriptomics

Differentially expressed genes across each clone were analyzed using the Seurat R package (version 4.0.5). Space Ranger output files were imported, and the data were normalized and scaled using the default NormalizeData and ScaleData functions. Differential gene expression analyses were performed comparing groups using the FindMarkers function, with the following parameters: log (fold change) ≥ 0.1 and adjusted P-value ≤ 0.05 (Fig. 4G). The list of CHIP-seq-based STAT1 target genes was downloaded from previous literature 102. For differentially expressed genes in the CRC4 slice, Gene Ontology (GO) enrichment analysis was performed by clusterProfiler (v4.6.0) (Fig. 4G) 103.

TCGA-COAD data analysis

Gene expression and clinical data for Colon Adenocarcinoma (COAD) were downloaded from The Cancer Genome Atlas (TCGA). Average copy number data were downloaded from the UCSC Xena browser (https://xena.ucsc.edu). For CNV analysis, only the samples with CNV information were obtained. The average log FC_CNV > 0.35 was defined as copy number gain, while the average log FC_CNV < -0.35 was defined as copy number loss (Fig. 4F). The signatures of each cell subpopulation were obtained in CRC single-cell data by the FindAllMarkers function in the Seurat package with default parameters. The normalized expression of the top 10 signatures in each cell subpopulation was used to calculate the average signature score. The reads per million (RPM) of known human endogenous retroviruses (HERV) in TCGA-COAD were downloaded from a previous report 104, 105. For HERV analysis, samples without HERV expression information were excluded. Overall HERV-H expression was calculated as the average expression of known transcribed HERV-H loci (HERV-H1-7) (Fig. 4K) 104. For chromosomal instability assessment, samples were stratified into tertiles by aneuploidy score, followed by comparative analysis of T cell marker expression across groups. The pathway activity of each sample was calculated using GSVA with default parameter settings (fig. S5K) 106. Survival analyses were then performed with the survival R package (version 3.5-3) (Fig. 4L and fig. S5R).

MSK-IMPACT data analysis

The genomic and clinical data of 71 CRC patients with Immune Checkpoint Blockade (ICB) therapy were downloaded from published studies 51, 52. The copy number status of each loci was determined using the GISTIC v2.0 function in the GenePattern platform 107, 108. The average log FC_CNV > 0.35 was defined as copy number gain, while the average log (copy number/2) < -0.35 was defined as copy number loss. Survival analyses were then performed with the survival R package (Fig. 4L and fig. S5, P to Q).

Experimental and Analytical Approach to Mouse Ear

Ear wound creation and treatment

All mouse experiments were conducted in accordance with all relevant ethical regulations and were approved by the Institutional Animal Care and Use Committees of Beijing Institute of Genomics, Chinese Academy of Sciences (approval number: 2019A010). MRL/MpJ (MRL; RRID: IMSR_JAX:000486) and C57BL/6 (B6; RRID: MGI:2159769) mice of both sexes, aged 8-10 weeks, were used for ear wound creation. All mice were housed under specific pathogen-free (SPF) conditions in the barrier facilities of the Tsinghua University Laboratory Animal Center and the Laboratory Animal Center of the Institute of Genetics and Developmental Biology, Chinese Academy of Sciences. Animals were maintained under a 12 h light/12 h dark cycle, at an ambient temperature of 22 ± 1 °C and a relative humidity of 50 ± 10%, with food and water provided ad libitum. For wound generation, mice were briefly anesthetized with isoflurane. The dorsal and ventral aspects of mouse ears were sanitized with cotton balls saturated in 75% alcohol. Mouse ears were delicately flattened using forceps, and a tissue-penetrating ear puncher (2 mm) was used to create perforating injuries within the vascularized regions of the ear. Hemostatic powder was applied to the wounds using disposable cotton swabs. All instruments were sterilized under high pressure before use.

To block the differentiation of Fmod+ fibroblasts, mice were either locally injected with rosiglitazone (3 mg/kg) every other day post-injury or smeared with rosiglitazone dissolved in petrolatum every day. For the injection administration, the control group received an equivalent volume of DMSO. For the topical administration, the control group received petrolatum containing DMSO.

To block IL11-IL11RA signaling, mice were locally injected with anti-IL11RA antibody (2.5 µg/wound) at day 3 post-injury. The control group received an equivalent volume of anti-rabbit IgG.

Flow cytometry

Single-cell suspensions were obtained from the ear tissues of B6 and MRL mice at day 7 post-injury. Specifically, mouse ear tissues were mechanically minced and then incubated at 37°C for 5-10 minutes with 40 µL collagenase II (50 mg/mL) and 40 µL collagenase IV (50 mg/mL). Following a brief low-speed centrifugation (300 × g) to precipitate undigested tissue fragments, the supernatant was transferred into a stop solution composed of fetal bovine serum, EDTA, and D-HBSS. Additional digestion enzymes were added to the partially digested tissue fragments, and the digestion process was repeated until complete tissue digestion was achieved. The digested mixture was filtered through a 70-µm membrane and centrifuged at 4°C at 1200 rpm for 10 minutes. After discarding the supernatant, the cell pellet was washed in FACS buffer (PBS supplemented with 2% FBS and 1 mM EDTA). The cells were first stained with Zombie Aqua™ Fixable Viability Kit (BioLegend Cat# 423102, 1:1000) to exclude dead cells, followed by incubation with Fc Block (BioXCell Cat# BE0307, clone 2.4G2, RRID:AB_2736987) at a 1 ug/mL concentration for 10 minutes. Subsequently, the following antibodies were used for surface staining: CD45-BUV395 (BD Biosciences Cat# 564279, clone 30-F11, RRID:AB_2651134), CD31-FITC (Invitrogen Cat# 11-0311-81, clone 390, RRID:AB_465012), CD29-PE (Invitrogen Cat# 12-0291-81, clone HMb1-1, RRID:AB_657733), Sca-1-APC/Cy7 (BioLegend Cat# 108125, clone D7, RRID:AB_10639725), CD34-Pacific Blue (BioLegend Cat# 119321, clone MEC14.7, RRID:AB_10900980), CD55-PE/Cy7 (BioLegend Cat# 131813, clone RIKO-3, RRID:AB_2800633), and CD90.2-BV711 (BioLegend Cat# 105349, clone 30-H12, RRID:AB_2800564). Antibodies were applied at a 0.2 ug/ml concentration. Staining was performed for 30 minutes on ice in the dark. Flow cytometry was performed using a BD Symphony A5 cytometer for analysis and a BD Fusion cell sorter for sorting experiments. The flow cytometry results were analyzed using FlowJo™ v10 Software (BDR Life Sciences).

Multi-color immunohistochemistry (mIHC)

At predetermined time points, murine ears were excised and fixed in 10% formalin for 48 hours at 4°C. The fixed tissues were then embedded in paraffin and sectioned at 7 µm thickness for histological evaluation. Tissue sections were then deparaffinized, rehydrated and stained using 7-color multispectral Opal reagents (AKOYA BIOSCIENCES, Cat# number NEL81 1001KT). Heat-mediated antigen retrieval was performed using a microwave for rehydrated tissue sections. Primary antibodies against CD11b (Abcam Cat# ab133357, clone EPR1344, RRID:AB_2650514; 1:500), IL-1β polyclonal antibody (Invitrogen Cat# P420B, RRID:AB_223478; 1:200), CD163 (Abcam Cat# ab182422, clone EPR19518, RRID:AB_2753196; 1:1000), CD68 (Proteintech Cat# 28058-1-AP, RRID:AB_2881049; 1:2000), CD74 (Abcam Cat# ab289885, clone EPR25399-94, 1:2000), GPNMB (Proteintech Cat# 66926-1-Ig, clone 2B10B8, RRID:AB_2882253; 1:200), α-SMA(Abcam Cat# ab124964, clone EPR5368, RRID:AB_11129103; 1:2000), PDGFR alpha (PDGFRα, Abcam Cat# ab96569, clone 16A1, RRID:AB_10687154; 1:500), CD34 (Invitrogen Cat# MA5-29674, clone 013, RRID:AB_2785500; 1:2500), OGN (Proteintech Cat# 12755-1-AP, polyclonal, RRID:AB_10646479; 1:200), SOX9 (CST Cat# 82630T, clone D8G8H, 1:500), Collagen Type II (Proteintech Cat# 28459-1-AP, polyclonal, RRID:AB_2881147; 1:200), IGFBP5 (Proteintech Cat# 55205-1-AP, polyclonal, RRID: AB_2736835; 1:200), PRG4 (Proteintech Cat# 25732-1-AP, polyclonal, RRID:AB_2880215; 1:200), K14 (Santa Cruz Cat# sc-53253, clone LL001, RRID:AB_2134820; 1:300), CD31 (Abcam Cat# ab281583, clone RM1006, RRID:AB_3096925; 1:2000), FcεRIα (Proteintech Cat#10980-1-AP, polyclonal, RRID: AB_2103033; 1:50), Ly6G (Abclonal Cat# A22270, RRID:AB_3718562; 1:200), IL11RA (Invitrogen Cat# PA5-47667, polyclonal, RRID:AB_2607330; 1:500) and Osteopontin (Proteintech Cat# 22952-1-AP, polyclonal, RRID: AB_2783651; 1:500) were individually applied. The secondary antibody Polymer HRP Ms + Rb (Cat# ARH1 001EA) was applied and incubated at room temperature for 10 minutes. Signal amplification was performed using the Opal Fluorophore Working Solution on the slides at room temperature for 10 minutes. Finally, nuclei were stained with DAPI (Invitrogen Cat# D1306, RRID: AB_2629482, 5 mg/mL, 1:2500). mIHC imaging was collected using the Vectra Polaris Imaging System (PerkinElmer) 100 at 20x magnification and analyzed using the InForm Advanced Image Analysis Software (PerkinElmer, version 2.4.1).

Cell segmentation and quantification were carried out using QuPath 101. In short, DAPI stains were used for cell detection with the following parameters: requestedPixelSize 0.5 μm, backgroundRadius 8 μm, sigma 0.5 μm, minArea 10 μm, maxArea 400 μm, threshold 0.5, cellExpansion 5 μm. Cell subpopulations were identified by training an object classifier based on the combination of the marker’s expression. In the epidermis, keratinocytes were measured, while in the dermis, a comprehensive identification and quantification of other cellular subpopulations were carried out. For the pseudo-spot visualization, we first divide the mIHC data into densely distributed bins, each with a width and length of approximately 50 μm. Then, we calculate the proportion of each cell subpopulation in the pseudo-spots by dividing the number of specified cell subpopulations by the total cell count defined by DAPI.

Single-cell RNA-seq

Single-cell suspensions were obtained from B6 and MRL mouse ear tissue at day 7 post-wounding. Calcein (Invitrogen Cat# C3099) and 7-AAD (Sigma-Aldrich Cat# 129935) were freshly added at final concentrations of 1 μM and 1 μg/mL, respectively, before cell sorting. Then FITC+ 7-AAD- single cells were sorted for library construction of scRNA-seq. The libraries were performed using Chromium Next GEM Single Cell 3’ GEM, Library & Gel Bead Kit V3 (PN-1000075) purchased from 10x Genomics according to the protocols provided by manufacturers. The libraries were sequenced on an Illumina Novaseq 6000 platform.

Spatial gene expression assay

Fresh mouse ear tissues were embedded in the Optimal Cutting Temperature Compound (OCT, Sakura Tissue-TEK), and stored at -80°C. The 10-μm section was placed on the pre-chilled Optimization slides (Visium, 10x Genomics, PN-1000193). After collecting histological morphological information, the tissue was processed within the optimal permeabilization time of 18 minutes to digest and release mRNA from tissue slides. The released mRNA was subsequently captured by proximal oligonucleotides within the tissue, culminating in the generation of the respective libraries. The libraries were sequenced on an Illumina HiSeq X Ten platform.

Single-cell RNA-seq data processing

10x Genomics scRNA-seq data were aligned and quantified using the Cell Ranger Single-Cell Software Suit (version 3.0.2) against the mm10 reference genome. Genes with fewer detected cell numbers (fewer than 3) and cells with fewer detected genes (fewer than 200) were removed. To further filter out low-quality data, the following criteria were applied to filter cells: (1) the number of expressed genes greater than 500 or less than 5,000; (2) mitochondrial proportion less than 20%. To cluster single cells by their expression, the unsupervised graph-based clustering algorithm implemented in the Seurat R package (version 3.1.4) was used. The top 2,000 variable genes identified by the vst method were used to perform principal component analysis. The RunUMAP and RunTSNE functions were employed to generate UMAP and t-SNE dimensionality reduction results, while the FindNeighbors function was used to construct a shared nearest neighbor graph. The DoubletFinder was utilized to predict potential doublets and filter them out.

Spatial transcriptomics data processing

Using the Space Ranger software (https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/tutorials), the 10x Visium spatial transcriptomics data was aligned and quantified against the mm10 version of the mouse reference genome. Downstream analysis was conducted using the R package Seurat. Integration of unpaired scRNA-seq data with spatial transcriptomics sequencing data was performed for deconvolution analysis to obtain the proportion of each cell subtype within each spot of the spatial transcriptome. Parameters for UCASpatial: downsample_n = 1000, meta.resolution = 200, meta.filter = 0.95, cos.filter.threshold = 0.1, the others use the default. Parameters for other methods: Default. Using the threshold of 3% to filter the noise after decomposition by each method (fig. S6C).

Spatial transcriptomics data visualization

The original data placed mouse ear tissues on two slides: tissues from day 0, day 7, and day 15 were on the first slide, and tissues from day 3 were on the second slide. To better visualize the spatial transcriptomics data of mouse ears, we merged all the tissues into one slide and aligned them based on the time point post-wounding. Specifically, we first extracted the coordinates of each spot in all 4 tissues. We then transformed the coordinates of spots from day 3, with each row replaced by the column, and each column replaced by the row. We next transformed the coordinates of each spot’s columns, abiding by the time points. ggplot2 was used for visualizing spatial transcriptomics data after coordinate transformation.

The murine ear tissues have been annotated into four groups of regions based on H&E staining: unwounded (day 0), wound bed (1 mm radius), wound proximal (0-1 mm from the wound bed), and wound distal (> 1 mm from the wound bed).

Functional enrichment analysis in murine wound healing data

The signatures of macrophage subpopulations, Fmod+ fibroblast, and Igfbp5+ chondrocyte in the murine wound healing dataset were obtained in scRNA-seq data by the FindAllMarkers function in the Seurat package with min.pct = 0.1, min.diff.pct = 0.1, logfc.threshold = 0.25 and other default parameters. Metascape was used to perform functional enrichment analysis for Fmod+ fibroblast and Igfbp5+ chondrocyte (Fig. S7, L and M) 109. The enrichKEGG in the clusterProfiler package was used to perform functional enrichment analysis for macrophage subpopulations (fig. S6N) 103.

Trajectory interference for chondrocyte subpopulations

Monocle3 was used to infer the trajectory of chondrocytes in scRNA-seq data from the murine wound healing model (fig. S7I) 110-113. The top 20 PCA defined in scRNA-seq of chondrocytes were used. The reduce_dimension function was used to reduce the dimensions of chondrocytes (Parameters: umap.metric = “cosine”, umap.min_dist = 0.05, umap.n_neighbors = 20). cluster_cells function with parameters cluster_method = “Louvain”, resolution = 0.001, random_seed = 123 was used to cluster the cells. The learn_graph function with default parameters was used to learn the principal graph from the reduced dimension space using reversed graph embedding.

Ligand and receptor interaction analysis

To identify the significant and specific ligand-receptor pairs between Igfbp5+ chondrocyte, Fmod+ fibroblast, and Cd36+ Gpnmb+ Il1b- macrophage, the CellPhoneDB tool was used in scRNA-seq data from murine wound healing 114, 115. Cell expression was normalized by dividing the total counts. The interaction pairs with P-value < 0.01 or the mean attraction strength > 0 were retained.

Statistics & Reproducibility

No statistical method was used to predetermine sample size. No data were excluded from the analyses. Mice were randomly assigned to experimental groups. The Investigators were not blinded to allocation during experiments and outcome assessment. Sample size and statistical tests are described in the figure legends. Data are expressed as mean ± SEM unless otherwise specified. Statistical tests were performed using R statistical software.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The current manuscript presents a strong case for the performance of UCASpatial relative to other algorithms designed for the deconvolution of spot-based spatial transcriptomics platforms. The authors have conducted benchmarking and validation using orthogonal techniques to substantiate their findings and have provided a comprehensive response to the previous round of review. Overall, the methodology appears robust following the earlier revisions. That said, I would encourage the authors to moderate some of their claims, particularly given that the colorectal cancer (CRC) data were derived from a very limited number of samples.\nSpecific comments:\n1. In Figure 3F, it appears that the authors suggest they can identify Tem, Tex, and Trm cells using PD-1, CD103, and GZMK. However, these markers alone are insufficient for such precise identification. More generally, the correspondence between cell populations defined by scRNA-seq and the surrogate markers employed in multiplex IHC should be interpreted and presented with greater caution.\n2. In the discussion of simulated data, the authors state that simulated spots reflect the actual cellular composition found in “typical tumor architectures.” It is unclear what is meant by “typical” in this context, as tumor architecture is highly variable. The only generalizable pattern might be the inverse relationship between epithelial and stromal cell frequencies across tumor core and stromal regions, respectively.\n3. Some sections require linguistic revision, as they currently read awkwardly. For instance, the repeated use of the expression “the colorectal cancer” should be corrected. Similarly, terms such as “infiltrating tumor clones” are ambiguous—do the authors mean tumor regions composed of different clonal populations? Clarification would improve readability and precision.", "reviewer_2": "The manuscript presents UCASpatial, an entropy-weighted regression framework for deconvolving spot-based spatial transcriptomics and demonstrates its utility in human colorectal cancer and in murine ear-wound regeneration models. The study is carefully executed and clearly written. On the wound-healing side, the authors provide extensive spatial and multiplexed immunohistochemistry (mIHC) validation across time points and complement the deconvolution with in vivo perturbations (rosiglitazone and anti-IL11RA). These data convincingly support the central claims that (i) UCASpatial resolves low-abundance, transcriptionally subtle immune subsets and (ii) a chondrocyte-macrophage–fibroblast triad, engaging IL11–IL11RA signalling, contributes to scarring in C57BL/6 mice.\nBased on the current version, the authors appear to have addressed prior reviewer concerns with substantial additional data and clearer benchmarking/validation. I do not see unresolved major issues that would preclude publication. I recommend a minor revision, focused on small corrections noted below.\nPage 12, lines 492–493: the text refers to the regeneration stage, so Fig. 5H (regeneration) should be cited instead of Fig. 5G (return-to-unwounded).\nPage 12, lines 502–504: the manuscript states: “We identified a scarring-associated myofibroblast subpopulation, Lrrc15⁺ fibroblasts, located within the wound bed at the return-to-unwounded stage in both MRL and B6 mice (Fig. S6A–B).” Please support the “scarring-associated” designation by citing prior work linking LRRC15⁺ myofibroblasts to fibrosis/scar formation.\nOverall, the manuscript is strong and, with the small corrections and clarifications above, will make a valuable contribution to spatial transcriptomics and regenerative biology.", "reviewer_3": "The authors have addressed most of the concerns raised by previous reviewers and have made improvements to the manuscript. Major additions include Python integration, comprehensive experimental validation, and expanded performance benchmarking. However, several concerns require further clarification to strengthen the methodological foundation and ensure broad applicability across different platforms and biological contexts:\n(Reviewer #3) Is there any benchmarking on the meta-purity filtering with respect to whether it leads to the loss of certain cell states because of the filtering on intermediate cell subsets? Loss of these cell types might bias biological conclusions. The authors need to be explicit about what is getting filtered out and the impact.\nThe reviewer specifically asked about \"certain cell states\", but the authors only evaluated predefined, discrete cell types in their analysis. A PCC >0.9 for overall expression profiles does not prove preservation of rare cell states or cellular heterogeneity within cell types, as this metric averages across all cells and could easily mask the loss of small populations of transitional cells. Furthermore, the authors provide no analysis of what specific cells or states are filtered out at their arbitrary 95% threshold, no biological justification for why this threshold was chosen, and no sensitivity analysis showing how different thresholds affect both computational performance and biological interpretation. In addition, they did not evaluate whether the filtered cells represent biologically meaningful intermediate states.\n(Reviewer #3) How would biases due to dissociation (e.g., loss of rare or fragile cell types) affect the cGEP-based deconvolution? Can it distinguish their absence?\nThe author conducted simulations removing only CD8+ Tex and fibroblasts, showing maintained performance but acknowledged that \"UCASpatial struggles to detect the loss of rare or fragile cell types in scRNA-seq reference.\" This response only tested two specific cell types rather than systematically evaluating how reference quality affects performance across different tissue types and biological contexts. In addition, the author should improve the method by providing guidelines for users to assess reference data quality and detect when reference data is inadequate for their specific tissue context.\n(Reviewer #4) A recently published ST deconvolution tool, Spotiphy (Nature Methods 2025, Jiyuan Yang &Ziqian Zheng, PMID: 40074951), proposed a strategy for selecting informative genes for each cell type from scRNA-seq data, UCSpatial is conceptually similar to Spotiphy in terms of cGEP identification. A more thorough comparison of the performance between Spotiphy and UCSpatial is warranted.\nThe reviewer requested a more thorough comparison with Spotiphy. The authors added Spotiphy to some benchmarking analyses but provided only superficial comparisons limited to basic performance metrics. The comparison lacks detailed analysis of methodological differences between the entropy-based weighting approach in UCASpatial versus Spotiphy's informative gene selection strategy. In addition, the author can address the concern of reviewer #4 by discussing scenarios where each method might be preferred over the other.\n(Reviewer #4) In the benchmarking experiments, most evaluation metrics lack direct ground-truth support. The Slide-seq V2 dataset, which includes ground-truth labels, was not used for benchmarking UCSpatial against other deconvolution tools (only compared to RCTD) (Figures S1AB-AD). The authors need to include a comparative analysis between the predictions from these algorithms and the ground-truth labels in the Slide-seq V2 mouse hippocampus dataset.\nThe authors provided mIHC validation for immune tissues but failed to utilize the available Slide-seq V2 ground-truth data for systematic benchmarking. While they justify focusing on immune microenvironments, this does not excuse missing the opportunity to leverage established hippocampus cell type annotations for objective method comparison. The authors briefly mentioned testing on Visium HD and Slide-seq V2 data but provided minimal benchmarking results (only RCTD was compared). In addition, the author only showed spatial distribution of cell subpopulations obtained from UCASpatial, while the spatial plots obtained from RCTD were missed.\nMajor Concerns in the Computational Methodology Section\nWhile the overall methodology is well-implemented, several aspects could benefit from additional theoretical justification and expanded validation to enhance the robustness and generalizability of the approach:\n1. Line 81 and line 718, Shannon entropy: The Shannon-entropy-based weighting strategy lacks clear justification for why Shannon entropy is the best metric for determining gene importance. While the authors demonstrated UCASpatial's performance on in real murine wound healing and human CRC dataset to address Reviewer #1's high noise concerns, they should compare their approach against simpler alternatives like variance-based or mutual information weighting to show that Shannon entropy provides the best results. Without such comparisons, the decision to use Shannon entropy in UCASpatial seems somewhat arbitrary and weakens the methodological foundation.\n2. Line 134, Comparison Scope: The comparison with recently published methods like Spotiphy is limited to basic performance metrics. A more comprehensive analysis of why (by utilizing what information that alternatives ignore) UCASpatial outperforms alternatives is needed.\n3. Line 791 and Table S1, Simulation strategy: Despite adding mIHC validation, the benchmarking still relies heavily on simulated data. However, the simulation strategy may not capture the full complexity of real spatial transcriptomics data, such as the spatial correlation structures between neighboring spots and cell-type-specific capture efficiency variations. The author should consider these features to construct more reliable synthetic spots."}, "author_response": "Point-by-point response\nReviewers' Comments:\nReviewer #1:\nThe current manuscript presents a strong case for the performance of UCASpatial relative to other\nalgorithms designed for the deconvolution of spot-based spatial transcriptomics platforms. The\nauthors have conducted benchmarking and validation using orthogonal techniques to substantiate\ntheir findings and have provided a comprehensive response to the previous round of review. Overall,\nthe methodology appears robust following the earlier revisions. That said, I would encourage the\nauthors to moderate some of their claims, particularly given that the colorectal cancer (CRC) data\nwere derived from a very limited number of samples.\nSpecific comments:\n1.In Figure 3F, it appears that the authors suggest they can identify Tem, Tex, and Trm cells using PD-\n1, CD103, and GZMK. However, these markers alone are insufficient for such precise identification.\nMore generally, the correspondence between cell populations defined by scRNA-seq and the surrogate\nmarkers employed in multiplex IHC should be interpreted and presented with greater caution.\nResponse: We thank the reviewer for highlighting the need for precision and caution in correlating\nscRNA-seq clusters with protein markers. We agree that the correspondence between scRNA-seq-\ndefined clusters and protein markers used in mIHC is not one-to-one and requires careful interpretation.\nTo address this concern, we have revised the manuscript to adopt a more conservative and descriptive\nnomenclature for the mIHC-defined CD8+ T cell subpopulations. Specifically, we now refer to them\nas \"PD-1+ CD8+ T cells,\" \"CD103+ CD8+ T cells,\" and \"GZMK+ CD8+ T cells,\" rather than definitively\nlabeling them as CD8+ Tex, Trm, and Tem cells, respectively (Pages 7-8, lines 303-307; Fig. 3F-G).\nTo ensure consistency, we have also updated the figure legends for all relevant supplementary figures\n(e.g., fig. S4C and fig. S6D) to explicitly state that these marker combinations are used as surrogates\nto represent their corresponding cell states. Consequently, we have changed the wording in our\nconclusions from \"validated\" to \"supported\" to more accurately reflect that the mIHC data provides\ncorrelative and supportive evidence for our findings. We believe these changes present our findings\nwith the appropriate level of caution.\nFigure 1 for Reviewer #1. Nomenclature for the mIHC-defined CD8+ T cell subpopulations.\n(A) Left, representative immunofluorescence imaging of CRC tissue (CRC9). The staining images indicate\nT cell subpopulations within CRC9. The cyan, green, and magenta colors indicate PD-1+ CD8+ T cells, CD103+\nCD8+ T cells, and GZMK+ CD8+ T cells, respectively. The dashed white lines indicate the tumor-invasive front.\nScale bar, 100 μm. (B) Left, boxplot showing the distance between CD8+ T cell subpopulations and tumor\nregions in 10 ROIs. Right, bar plot showing the proportion of CD8+ T cell subpopulations within the tumor\nregion in 10 ROIs. ROI, region of interest.\n2. In the discussion of simulated data, the authors state that simulated spots reflect the actual cellular\ncomposition found in “typical tumor architectures.” It is unclear what is meant by “typical” in this\ncontext, as tumor architecture is highly variable. The only generalizable pattern might be the inverse\nrelationship between epithelial and stromal cell frequencies across tumor core and stromal regions,\nrespectively.\nResponse: We thank the reviewer for this valuable comment. We recognize that the phrase “typical\ntumor architectures” was imprecise, given the high heterogeneity of tumor spatial organization and the\ndifficulty of defining a single standard pattern. To address this, we have revised the relevant statement\nto: “The cellular composition of simulated spots was designed to recapitulate the spatially correlated\ncell-type distribution patterns observed in tumor microenvironments, particularly the relationship\nbetween epithelial, stromal, and immune cell frequencies across tumor core, margin, and stromal\nregions.” This wording more accurately reflects the rationale behind our simulation design and avoids\nthe potential misinterpretation associated with the term “typical.” We appreciate the reviewer’s\ninsightful suggestion and have updated the manuscript accordingly.\n3. Some sections require linguistic revision, as they currently read awkwardly. For instance, the\nrepeated use of the expression “the colorectal cancer” should be corrected. Similarly, terms such as\n“infiltrating tumor clones” are ambiguous—do the authors mean tumor regions composed of different\nclonal populations? Clarification would improve readability and precision.\nResponse: We thank the reviewer for their attentive reading and for the valuable suggestions to\nimprove the linguistic quality and precision of our manuscript. We have undertaken a thorough revision\nin direct response to the points raised. Specifically, the repeated expression “the colorectal cancer” has\nbeen corrected by removing the definite article for more idiomatic phrasing. The ambiguous phrase\n‘infiltrating tumor clones’ has been clarified to accurately reflect our intended meaning, now reading\n‘immune cells infiltrating into distinct tumor clone regions’. In addition, we have performed a\ncomprehensive check of the manuscript for grammar, spelling, and punctuation. Several such errors\nwere identified and corrected (e.g., Lines 385, 667-669). We believe these revisions have significantly\nenhanced the clarity and readability of the manuscript.\nReviewer #2:\nThe manuscript presents UCASpatial, an entropy-weighted regression framework for deconvolving\nspot-based spatial transcriptomics and demonstrates its utility in human colorectal cancer and in\nmurine ear-wound regeneration models. The study is carefully executed and clearly written. On the\nwound-healing side, the authors provide extensive spatial and multiplexed immunohistochemistry\n(mIHC) validation across time points and complement the deconvolution with in vivo perturbations\n2\n(rosiglitazone and anti-IL11RA). These data convincingly support the central claims that (i)\nUCASpatial resolves low-abundance, transcriptionally subtle immune subsets and (ii) a chondrocyte-\nmacrophage–fibroblast triad, engaging IL11–IL11RA signalling, contributes to scarring in C57BL/6\nmice.\nBased on the current version, the authors appear to have addressed prior reviewer concerns with\nsubstantial additional data and clearer benchmarking/validation. I do not see unresolved major issues\nthat would preclude publication. I recommend a minor revision, focused on small corrections noted\nbelow.\n1. Page 12, lines 492–493: the text refers to the regeneration stage, so Fig. 5H (regeneration) should\nbe cited instead of Fig. 5G (return-to-unwounded).\nResponse: We thank the reviewer for carefully noting this citation error. In the revised manuscript, we\nhave corrected the figure reference on Page 12, line 514, to cite Fig. 5H, which corresponds to the\nregeneration stage, instead of the incorrect Fig. 5G (return-to-unwounded).\n2. Page 12, lines 502–504: the manuscript states: “We identified a scarring-associated myofibroblast\nsubpopulation, Lrrc15⁺ fibroblasts, located within the wound bed at the return-to-unwounded stage in\nboth MRL and B6 mice (Fig. S6A–B).” Please support the “scarring-associated” designation by citing\nprior work linking LRRC15⁺ myofibroblasts to fibrosis/scar formation.\nResponse: We thank the reviewer for this helpful comment. To support our designation of the Lrrc15⁺\nmyofibroblast population as scarring-associated, we now emphasize previously published evidence\nrather than referencing our own figures. Prior studies have consistently shown that Lrrc15⁺ fibroblasts\nrepresent a fibrotic or scar-associated myofibroblast state across diverse pathological contexts,\nincluding idiopathic pulmonary fibrosis and ulcerative colitis1-3. These works also demonstrate that\nLrrc15⁺ fibroblasts are enriched for canonical ECM-remodeling and scar-forming gene programs (e.g.,\nPostn, Thy1, Acta2, Col1a1, Col1a2), supporting their functional alignment with scarring phenotypes1-\n3.\nIn the revised manuscript, we have added citations to these supporting studies (Page 12, line 525;\nreferences 73–75), including reports by Franzen et al., Abedini et al., and Buechler et al.\nOverall, the manuscript is strong and, with the small corrections and clarifications above, will make\na valuable contribution to spatial transcriptomics and regenerative biology.\nResponse: We sincerely thank the reviewer for the encouraging and supportive assessment of our\nmanuscript. We greatly appreciate the thoughtful suggestions and comments throughout the review\nprocess, which have helped us improve the clarity and quality of our work. We hope that our revisions\neffectively address the reviewer’s concerns and that our study will make a meaningful contribution to\nthe fields of spatial transcriptomics and regenerative biology.\nReviewer #3:\nThe authors have addressed most of the concerns raised by previous reviewers and have made\nimprovements to the manuscript. Major additions include Python integration, comprehensive\n3\nexperimental validation, and expanded performance benchmarking. However, several concerns\nrequire further clarification to strengthen the methodological foundation and ensure broad\napplicability across different platforms and biological contexts:\n1. ‘(Reviewer #3) Is there any benchmarking on the meta-purity filtering with respect to whether it\nleads to the loss of certain cell states because of the filtering on intermediate cell subsets? Loss of\nthese cell types might bias biological conclusions. The authors need to be explicit about what is getting\nfiltered out and the impact.’\nThe reviewer specifically asked about \"certain cell states\", but the authors only evaluated predefined,\ndiscrete cell types in their analysis. A PCC >0.9 for overall expression profiles does not prove\npreservation of rare cell states or cellular heterogeneity within cell types, as this metric averages\nacross all cells and could easily mask the loss of small populations of transitional cells. Furthermore,\nthe authors provide no analysis of what specific cells or states are filtered out at their arbitrary 95%\nthreshold, no biological justification for why this threshold was chosen, and no sensitivity analysis\nshowing how different thresholds affect both computational performance and biological interpretation.\nIn addition, they did not evaluate whether the filtered cells represent biologically meaningful\nintermediate states.\nResponse: We thank the reviewer for this important and insightful comment. We acknowledge that\ndefining “transitory” or intermediate cell states is inherently challenging, as these cells often exhibit\ncontinuous transcriptional gradients rather than discrete boundaries. To address this, we designed the\nmeta-purity filtering strategy aiming to prioritize the robustly extract features of well-defined cell types.\nThis approach aggregates single cells into “meta-cells” and filters out those with mixed composition\n(i.e., purity below a given threshold), thereby reducing the impact of ambiguous signals on the\nidentification of cell-identity-specific gene expression profiles.\nRegarding the default purity threshold of 0.95, this value applies specifically to the purity of meta-\ncells, not directly to individual cell types. Building on the preliminary benchmarks in our original\nsubmission (fig. S1Q-S, now fig. S1S-U in revised manuscript), we have conducted systematic\nbenchmark analysis using simulated datasets to assess the effects of varying this hyperparameter on\ndeconvolution performance (Figures 1A-C for Reviewer #3). Through these analyses, the 0.95\nthreshold was found to provide relatively stable and optimal results, balancing filtering stringency and\nretention of meaningful biological signals.\nImportantly, we implemented meta-purity filtering as an optional, user-controlled parameter.\nUsers can adjust the threshold or disable this filter to preserve maximal cellular heterogeneity,\ndepending on their specific biological questions.\nWe have included these clarifications and relevant benchmarking results in the revised manuscript\non fig. S1S-U and the Meta-purity filtering in Methods, to provide transparency about the rationale\nand flexibility of our meta-purity filtering process.\n(A-C) Dot plots displaying the performance (A) PCC, (B) RMSE, and (C) F1 score with different\nhyperparameter values testing within the simulated ST dataset. Left panel, the performance of UCASpatial with\nentropy-based weight. Right panel, the performance of UCASpatial without entropy-based weight. Columns in\neach panel denote performance when varying selected hyperparameters (resolution and purity of meta-purity\nfiltering process) of UCASpatial. Red arrows denote the default hyperparameter values used in the benchmark\npart (resolution = 100, purity = 0.95). The red dashed line represents the performance (F1 score, PCC, and\nRMSE) of UCASpatial without the meta-purity filter.\n2. (Reviewer #3) How would biases due to dissociation (e.g., loss of rare or fragile cell types) affect\nthe cGEP-based deconvolution? Can it distinguish their absence?\nThe author conducted simulations removing only CD8+ Tex and fibroblasts, showing maintained\nperformance but acknowledged that \"UCASpatial struggles to detect the loss of rare or fragile cell\ntypes in scRNA-seq reference.\" This response only tested two specific cell types rather than\nsystematically evaluating how reference quality affects performance across different tissue types and\nbiological contexts. In addition, the author should improve the method by providing guidelines for\nusers to assess reference data quality and detect when reference data is inadequate for their specific\ntissue context.\nResponse: We thank the reviewer for raising this critical challenge. We fully agree that biases\nintroduced during scRNA-seq dissociation, such as the complete omission of rare or fragile cell types,\nconstitute an open and fundamental challenge for all deconvolution or mapping methods relying on\nreferences in this field. Our prior analysis was intended to illustrate the concept but was indeed limited\nin scope.\nguidelines, we have significantly enhanced our methodology by developing a Reference Diagnosis\nTool. This new feature allows users to assess the quality and sufficiency of their scRNA-seq reference\nprior to and during deconvolution (Figures 2A-C for Reviewer #3).\nSpecifically, we generated a visualization interface that helps evaluate whether the user-provided\nreference contains unique features for each cell subset under our feature extraction framework (Figure\n2 for Reviewer #3). We use a dot plot to visualize the correspondence between extracted topics and\ncell subset features, where an ideal scenario shows a diagonal distribution, indicating one-to-one\nmatching (Figure 2A for Reviewer #3). For example, in our scRNA-seq data used in the benchmark\npart, the algorithm effectively captures the core expression patterns of each cell subset, which is\nreflected in the clear diagonal correspondence in the dot plot (Figure 2A for Reviewer #3). If a topic\ncorresponds to multiple subpopulations, this may signal that the subset definitions need reevaluation\n(Figure 2B for Reviewer #3). For instance, we generated a “low-quality” reference by applying random\nperturbations (random seed = 123) to the original labels of the three CD8⁺ T-cell subpopulations—\n‘CD8_Tem’, ‘CD8_Tex’, and ‘CD8_Trm’— while all other cell-type labels remained untouched. After\nre-running the same mapping and visualization pipeline, the resulting dot plot no longer exhibited the\noriginal one-to-one correspondence; instead, multiple topics showed comparable, intermediate weights\nacross subpopulations, resulting in the loss of the one-to-one matching pattern in the dot plot (Figure\n2B for Reviewer #3). This demonstrates how inaccuracies or ambiguities in cell subset annotation can\nbe detected using this visualization tool.\nAdditionally, we recommend that when users lack suitable reference datasets for their tissue\ncontext, they consider using well-defined, publicly available high-quality scRNA-seq data. Our\nmethod supports the use of non-paired high-quality references, expanding applicability.\nThis tool provides a practical and systematic guideline for users: if the diagnosis indicates a poor\nmapping, the reference requires refinement (e.g., re-clustering, removing low-quality cells, or\nacquiring an alternative reference), thereby fulfilling the reviewer's request for a robust mechanism to\ndetect when reference data is inadequate for a specific tissue context. We greatly appreciate the\nreviewer’s insightful feedback, and we have integrated this feature into the UCASpatial package and\nfully documented it in the Methods section and the user tutorial.\nFigure 2 for Reviewer #3. Reference-quality diagnosis via topic-subset correspondence.\n(A) Ideal reference scenario: dot-plot matrix of 20 topics (y-axis) versus 20 cell subpopulations (x-axis) produced\nfrom the high-quality scRNA-seq reference. Dot size shows the weight magnitude of a topic in a cell subset, while\ncolor intensity indicates the exact weight value. A sharp diagonal pattern demonstrates one-to-one matching,\nconfirming that the latent topics uniquely capture the transcriptional fingerprints of individual cell subpopulations.\n6\n(B) Low-quality reference scenario: identical pipeline applied after randomly permuting the labels of the three CD8⁺\nT-cell subpopulations (CD8_Tem, CD8_Tex, CD8_Trm; random seed = 123). The diagonal structure collapses;\nmultiple topics now exhibit comparable weights across subpopulations, providing an immediate visual alert that\nsubpopulation definitions are ambiguous or mis-annotated.\n3. (Reviewer #4) A recently published ST deconvolution tool, Spotiphy (Nature Methods 2025,\nJiyuan Yang &Ziqian Zheng, PMID: 40074951), proposed a strategy for selecting informative genes\nfor each cell type from scRNA-seq data, UCSpatial is conceptually similar to Spotiphy in terms of\ncGEP identification. A more thorough comparison of the performance between Spotiphy and\nUCSpatial is warranted.\nThe reviewer requested a more thorough comparison with Spotiphy. The authors added Spotiphy to\nsome benchmarking analyses but provided only superficial comparisons limited to basic performance\nmetrics. The comparison lacks detailed analysis of methodological differences between the entropy-\nbased weighting approach in UCASpatial versus Spotiphy's informative gene selection strategy. In\naddition, the author can address the concern of reviewer #4 by discussing scenarios where each\nmethod might be preferred over the other.\nResponse: We thank the reviewer for highlighting the importance of a thorough comparison between\nUCASpatial and Spotiphy. We agree that positioning our tool against such state-of-the-art methods is\ncritical for demonstrating its utility.\nRegarding the performance benchmarking: We would like to clarify that comprehensive\nbenchmarking against Spotiphy was actually integral to our original submission. We have included\nSpotiphy in multiple evaluations across both simulated datasets (original Fig. 2A-G and fig. S1) and\nreal biology contexts (original fig. S3I-L, fig. S4E-H, and fig. S6G-I). We apologize if these results\nwere not sufficiently prominent or easy to locate in the previous version. To address this, we have now\nconsolidated these comparisons and made them more explicit in Figures 3-5 for reviewer #3.\nThese evaluations assessed multiple performance aspects beyond basic metrics (Figure 3 for\nReviewer #3), such as evaluation of cell-subpopulation-level PCC in human colorectal cancer and\nmurine wound healing (Figures 4A, B, E, and F for Reviewer #3), spatial pattern reconstruction of\nTLS in human colorectal cancer (Figures 4C and D for Reviewer #3). Consistently in both in-silico\nand in-real, we noticed that UCASpatial still exhibited the best performance, while Spotiphy exhibited\na 5/8 rank in our extensive benchmarking, showing a relatively poor performance in ST deconvolution\nat both spot-level and cell-type-level (Figures 3-4 for Reviewer #3). Furthermore, in the mIHC-based\nvalidation, UCASpatial showed markedly improved concordance with experimentally measured\nspatial cell-type densities. Specifically, UCASpatial more achieves more precise annotation of 12\nimmune cell subpopulations than existing methods in human colorectal cancer (Figures 5A-D for\nReviewer #3), as well as the dynamic redistribution of neutrophils, monocytes, and macrophage\nsubpopulations across the inflammatory-to-resolution phases in the wound-healing model (Figures 5E-\nF for Reviewer #3). These results indicate that UCASpatial not only performs well on simulated and\nreal ST datasets but also aligns more faithfully with high-resolution imaging ground truths, particularly\nin biologically complex contexts such as TLS formation and myeloid-state transitions during tissue\nrepair. We apologize if these analyses were not sufficiently prominent, and in the revised manuscript\nwe have made these comparisons more explicit and easier to locate.\nRegarding the methodological differences (Feature Selection vs. Entropy Weighting): We\n7\nappreciate the reviewer’s observation regarding the conceptual similarities in identifying informative\ngenes. However, there is a fundamental distinction for identifying informative genes. Spotiphy\n(feature-selection) adopt a “hard” selection strategy, retaining only the top informative genes and\ndiscards the rest. While this increases specificity, it may loss of subtle, low-abundance signals that are\ncritical for defining heterogeneous cell states. In contrast, UCASpatial (Entropy Weighting) adopts a\n\"soft\" weighting strategy. We do not discard genes; instead, we employ an entropy-based weighting\nscheme that adaptively highlights the most informative genes while still retaining the full gene space.\nThis enables UCASpatial to capture more nuanced and context-dependent cGEP variation without\nlosing potentially important signals. We have now clarified these methodological distinctions in the\nrevised manuscript.\nFinally, in response to the reviewer’s constructive suggestion, we have added a Discussion of\nscenarios in which each method may be preferred in our revised manuscript on Page15, lines 645-656:\n“Beyond UCASpatial’s strengths in resolving heterogeneous immune microenvironments and rare cell\nsubsets, other deconvolution tools retain niche advantages that justify their continued use. Spotiphy\nuniquely fuses histological images with expression to deliver single-cell-resolution whole-\ntranscriptome images. RCTD and cell2location perform comparably in benchmarks, yet RCTD’s\nlightweight Bayesian model demands few dependencies and minimal memory, making it attractive for\nrapid, resource-limited projects. Conversely, cell2location requires GPUs and large memory, but its\nseamless integration with Scanpy offers Python-native users a familiar, reproducible pipeline. CARD,\nby embedding neighborhood or morphological priors, is preferable for highly ordered tissues such as\ncortex or cardiac muscle where spatial continuity dominates. Thus, while UCASpatial maximizes\nperformance in heterogeneous, immune-rich microenvironments, method selection should balance\ncompute budget, programming ecosystem and the value of image information to match each biological\nquestion.”\nWe hope these clarifications and the expanded discussion fully address the reviewer’s concerns.\n(A)Graphic model showing the classification of 750 simulated spots into three groups, each representing a different\ncomplexity level - low, medium, and high, respectively. The classification is based on the entropy of cell-type\ncomposition for each spot, with each group containing 250 spots.\n(B-D) Boxplot showing spot-level root mean square error (RMSE) (B), PCC (C), and F1 score (D) for each method\nacross three complexity levels. Red boxes label the performance of Spotiphy. Boxplot elements are defined as: center\nline, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers. Statistical\nanalysis was performed using a two-tailed paired t-test. ns, non-significant, p > 0.05; *0.01 ≤ p < 0.05; **p < 0.01;\n***p < 0.001; ****p < 0.0001.\n(E) Line charts showing cell-subpopulation-level Pearson correlation coefficient (PCC) for each method across 20\ncell subpopulations. ‘avg’ and ‘maxm’ represent the average and maximum PCC, respectively, of the seven methods\nexcluding UCASpatial. The red arrow labels the performance of Spotiphy.\n(F)Left, line charts representing cell-subpopulation-level PCC of CD8+ effector memory T cells (CD8+ Tem), CD8+\nexhausted T cells (CD8+ Tex), CD8+ tissue-resident memory T cells (CD8+ Trm), and natural killer cells (NK) for\n9\neach method. Right, stacked bar plot showing the accumulated count of false predictions of these cell subpopulations\nfor each method. The red arrows label the performance of Spotiphy. Colors represent different cell subpopulations.\nFP, false positive; FN, false negative; FP + FN, total error frequency.\n(G) The cell-subpopulation-level PCC values for each method were binned into four confidence groups: 0-0.25 (low),\n0.25-0.5 (medium-low), 0.5-0.75 (medium-high), and 0.75-1 (high). Stacked bar plot showing the ratio of four\nconfidence groups for each method. The red arrow labels the performance of Spotiphy.\n(H) Graphic model representing the classification based on the proportion range of a specific cell subpopulation\nwithin the simulated spots. All simulated spots were classified into ten groups from proportion 0-10% to proportion\n0-100% for each cell subpopulation.\n(I)Line charts representing the average RMSE (left), PCC (middle), and F1 score (right) of 20 cell subpopulations\nacross various proportions for each method. Each dot in the line charts represents the average value from multiple\nsimulations. The red arrows label the performance of Spotiphy.\n(A) Boxplot showing cell-subpopulation-level PCC for 8 methods across T cell subpopulations, immune cell\nsubpopulations excluding T cells, non-immune subpopulations, and all cell subpopulations. Red boxes label the\nperformance of Spotiphy. Statistical analysis was performed using a two-tailed unpaired t-test. ns, non-significant,\np > 0.05; *0.01 ≤ p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. Boxplot elements are defined as: center line,\nmedian; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers.\n11\n(B) Heatmap showing cell-subpopulation-level PCC for 8 methods between estimated proportions and marker gene\nscores. UCASpatial_ablation represents UCASpatial without the entropy-based weight and meta-purity filter. The\nred arrow labels the performance of Spotiphy.\n(C) Stacked bar plots showing the composition of cell subpopulations in the TLS regions annotated by 8 methods.\nOnly slices containing TLS regions, as defined by pathological annotation, are displayed in the stacked bar plots. Red\nboxes label the performance of Spotiphy.\n(D) Spatial plot displaying the distribution of cell subpopulations within the TLS of the CRC5 slice, as annotated by\n8 methods. UCASpatial_ablation: UCASpatial without entropy-based weight and meta-purity filter. Red boxes label\nthe performance of Spotiphy.\n(E) Heatmap displaying cell-subpopulation-level Pearson correlation coefficient (PCC) between gene signature\nscores and estimated proportions from 8 methods. UCASpatial ablation represents the same framework as\nUCASpatial but with the core functions of entropy-based weighting and meta-purity filtering process disabled. The\nred arrow labels the performance of Spotiphy.\n(F) Boxplots showing cell-subpopulation-level PCC between gene signature scores and estimated proportions from\n8 methods across 23 cell subpopulations. Red boxes label the performance of Spotiphy. UCASpatial ablation\nrepresents the same framework as UCASpatial but with the core functions of entropy-based weighting and meta-\npurity filtering process disabled. The red dots represent the average PCC or rho of each method. Boxplot elements\nare defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points,\noutliers. Statistical analysis was performed using a two-tailed unpaired t-test. ns, non-significant, p > 0.05; *0.01 ≤ p\n< 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.\n(A) Left, bar plots displaying PCC and RMSE for each method in TLS regions for 12 immune cell subpopulations.\nRight, bar plots displaying PCC and RMSE for each method in tumor regions for 6 T cell subpopulations. Red boxes\nlabel the performance of Spotiphy.\n(B) Scatter plots showing the proportion of cell subpopulations evaluated by mIHC and each deconvolution method\nin tumor regions. The dashed line represents y = x. The proportions of cell subpopulations in tumor1 and tumor2\nwere calculated separately. Red boxes label the performance of Spotiphy.\n(C) Scatter plots showing the proportion of cell subpopulations evaluated by mIHC and each deconvolution method\nin TLS regions. The dashed line represents y = x. Red boxes label the performance of Spotiphy.\n(D) The proportion of cell subpopulations in the mIHC ROI and spatial transcriptomics evaluated by each\n13\ndeconvolution method. The mIHC data was separated into densely 50μm bins, and then the proportion of each\nsubpopulation was quantified (number of subpopulation cells / number of DAPI-identified cells). Completely gray\nsections represent that the proportion of indicated cell subpopulations in any spot is < 0.03. White, yellow, and red\ndashed lines indicate tumor1, tumor2, and TLS regions, respectively. Red boxes label the performance of Spotiphy.\n(E) Spatial distributions from ST data and representative immunofluorescence images of neutrophil (left), basophil\n(middle), and macrophage (right). The dashed rectangles indicate the wound bed. The dashed yellow lines indicate\nthe interface of the epidermis-dermis. The white dashed lines indicate the chondrocyte. Scale bar, 100 μm.\n(F) Left, bar plots showing the proportion of basophils in mIHC and ST data. For mIHC data, B6 day 3, n= 5;\nMRL day 3, n = 4; B6 day 7, n = 2; MRL day7, n = 1. Middle, bar plots showing the proportion of neutrophils in mIHC\nand ST data. For mIHC data, B6 day 0 and MRL day 0, n = 3; B6 day 3, n= 1; MRL day 3, n = 2; B6 day 7, n = 2;\nMRL day7, n = 1; B6 and MRL day 14, n = 3. Right, bar plots showing the proportion of macrophages in mIHC and\nST data. For mIHC data, B6 day 0 and MRL day 0, n = 1; B6 day 3, n= 1; MRL day 3, n = 2; B6 day 7, n = 3; MRL\nday7, n = 2. Red boxes label the performance of Spotiphy.\n4. (Reviewer #4) In the benchmarking experiments, most evaluation metrics lack direct ground-\ntruth support. The Slide-seq V2 dataset, which includes ground-truth labels, was not used for\nbenchmarking UCSpatial against other deconvolution tools (only compared to RCTD) (Figures S1AB-\nAD). The authors need to include a comparative analysis between the predictions from these\nalgorithms and the ground-truth labels in the Slide-seq V2 mouse hippocampus dataset.\nThe authors provided mIHC validation for immune tissues but failed to utilize the available Slide-seq\nV2 ground-truth data for systematic benchmarking. While they justify focusing on immune\nmicroenvironments, this does not excuse missing the opportunity to leverage established hippocampus\ncell type annotations for objective method comparison. The authors briefly mentioned testing on\nVisium HD and Slide-seq V2 data but provided minimal benchmarking results (only RCTD was\ncompared). In addition, the author only showed spatial distribution of cell subpopulations obtained\nfrom UCASpatial, while the spatial plots obtained from RCTD were missed.\nResponse: We sincerely thank the reviewer for highlighting the importance of this validation. We agree\nthat the Slide-seq V2 mouse hippocampus dataset serves as a critical gold standard for evaluating\ndeconvolution performance, particularly for assessing generalizability beyond immune-centric tissues.\n1. Following the reviewer’s suggestion, we have now performed a comprehensive benchmark on\nthe Slide-seq V2 dataset. Since strictly defined ground-truth cell type labels are not available\nexperimentally, we adopted the widely accepted evaluation strategy established by the CARD paper\n(Ma et al., Nat. Biotechnol., 2022)4. This involves verifying whether the deconvolved cell types\naccurately reconstruct known anatomical structures (e.g., CA1, CA3, Dentate Gyrus). Specifically, we\nvisualized the predicted top 1 cell type for each location (Figure 6A for Reviewer #3), and\nquantitatively evaluated the concordance between method-predicted cell-type proportions and known\nregion-specific marker expression in three anatomical structures4 by calculating Pearson correlation\ncoefficients (PCC). Our results show that UCASpatial achieves performance comparable to the best-\nperforming methods in terms of spatial resolution and accuracy (Figure 6B for Reviewer #3). We have\nalso provided the deconvolution results of all subpopulations from all methods (Figure 6C for\nReviewer #3).\nOur results show that UCASpatial is highly concordant with other methods, successfully captures\nkey anatomical structures, and confirms a comparable level of accuracy and spatial resolution on this\nneurobiological dataset. It is noteworthy that for CARD, we utilized a downsampled scRNA-seq\n14\nreference dataset consistent with the original CARD publication4 to faithfully reproduce its results,\nwhereas other methods were run using the full single-cell reference. Similar to the strategy in the\nCARD paper, we didn’t apply cell2location on this dataset due to its heavy computational burden. This\ndata has been added into our revised manuscript on Page 6, lines 220-229, and fig. S2F-H.\n15\n17\n(A) The top 1 cell type in each location inferred from seven different deconvolution methods. Colors represent\ndifferent cell types.\n(B) Bar plots showing the Pearson correlation coefficient between the marker gene expression and the predicted cell\ntype proportion region of each method. Wfs1, Cpne4, and C1ql2 represent CA1, CA3, and dentate gyrus, respectively.\n(C) Scatter plots showing the detailed deconvolution results of 17 cell types across spatial locations for seven methods.\nColor was shown to represent the minimal to maximal range of cell type proportions correspondingly.\n2. Regarding the reviewer's note on the limited benchmarking for Visium HD: We prioritized\nRCTD for comparison because other methods (e.g. Spotiphy, SPOTlight, stereoscope, CARD, and\ncell2location), faced compatibility issues or prohibitive computational demands with the 8µm high-\nresolution bins (which exceeds 500,000 total bins).\nAs requested by the reviewer, we have now included a direct side-by-side spatial heatmaps\ncomparison for all 34 cell subpopulations for both UCASpatial and RCTD (Figure 7 for Reviewer #3).\nWe acknowledge that this crucial data was not adequately highlighted in the original submission. We\nsincerely thank the reviewer for pointing this out, as it prompted us to find a better and more explicit\nway to display the superior performance of UCASpatial on high-resolution data. We confirm that these\nhigh-resolution images are still available on GitHub, but their inclusion in the revised manuscript is\nintended to ensure they receive the appropriate attention.\nFigure 7 for Reviewer #3. Benchmarking on the mouse hippocampus Slide-seq V2 dataset.\nSpatial plots displaying the distribution of all 34 cell subpopulations at spatial locations obtained from (A)\nUCASpatial and (B) RCTD in human CRC data based on 10X Visium HD platform with a bin size of 8μm × 8μm.\nColors indicate cell subpopulation proportions.\nWhile the overall methodology is well-implemented, several aspects could benefit from additional\ntheoretical justification and expanded validation to enhance the robustness and generalizability of the\napproach:\n1. Line 81 and line 718, Shannon entropy: The Shannon-entropy-based weighting strategy lacks clear\njustification for why Shannon entropy is the best metric for determining gene importance. While the\nauthors demonstrated UCASpatial's performance on in real murine wound healing and human CRC\ndataset to address Reviewer #1's high noise concerns, they should compare their approach against\nsimpler alternatives like variance-based or mutual information weighting to show that Shannon\nentropy provides the best results. Without such comparisons, the decision to use Shannon entropy in\nUCASpatial seems somewhat arbitrary and weakens the methodological foundation.\nResponse: We thank the reviewer for this thoughtful and constructive comment. We agree that\nproviding a clear justification for the choice of the weighting metric is essential for the methodological\nfoundation of UCASpatial.\nFirst, we chose Shannon entropy as a weighting metric because it offers specific advantages for\nscRNA-seq feature weighting compared to other metrics, such as mutual information (MI) and\nVariance-based methods. Shannon entropy quantitatively captures the specificity of gene expression\nacross cell subsets, providing good interpretability. Moreover, entropy calculation is computationally\nefficient and stable, making it well-suited for large-scale scRNA-seq data. In complex\nmicroenvironment contexts, many genes do not exhibit simple, exclusive high expression in a single\nsubpopulation; rather, they often show multimodal or subpopulation-specific partial expression\npatterns. Shannon entropy effectively reflects such expression variability and the evenness of\nexpression distribution. While mutual information (MI) can capture complex nonlinear dependencies\nand the association between gene expression and cell subpopulations, it typically requires large sample\nsizes and sophisticated estimation methods. In scRNA-seq data, where some subpopulations may\ncontain relatively few cells, MI estimation often suffers from instability and larger errors. Variance-\nbased weighting is straightforward and intuitive, highlighting expression variability, but is sensitive to\nextreme values and noise, and it cannot capture nonlinear dependencies.\nSecond, we conducted benchmarking experiments on simulated datasets comparing Shannon\nentropy weighting with both MI-based and variance-based weighting. Our results show that entropy\nweighting consistently outperforms the other two methods across multiple metrics, including Pearson\ncorrelation coefficient (PCC), root mean square error (RMSE), and F1 score (Figure 8 for Reviewer\n#3).\nBox plots displaying the PCC (left), RMSE (middle), and F1 score (right) of the deconvolution performance under\nthree different weighting strategies: entropy-, mutual information-, and variance-based, on the simulated datasets.\nBoxplot elements are defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5×\ninterquartile range; points, outliers. Statistical analysis was performed using a two-tailed paired t-test. **p < 0.01;\n***p < 0.001; ****p < 0.0001.\nAlthough Shannon entropy performed best in our benchmarks, we recognize that variance- and\nmutual information-based weighting schemes might be useful in specific user scenarios. Therefore,\nfollowing the reviewer’s suggestion, we have integrated these two alternative weighting options into\nUCASpatial, allowing users to select their preferred metric. We have updated the Methods section on\nPages 20-21, lines 822-842, to describe these new options.\n2. Line 134, Comparison Scope: The comparison with recently published methods like Spotiphy is\nlimited to basic performance metrics. A more comprehensive analysis of why (by utilizing what\ninformation that alternatives ignore) UCASpatial outperforms alternatives is needed.\nResponse: We thank the reviewer for the insightful comment. In addition to the expanded performance\ncomparison provided in our previous response (Figures 3-5 for Reviewer #3), we now include a clearer\nexplanation of why UCASpatial outperforms recently published alternatives such as Spotiphy in\nprinciple.\nExisting deconvolution methods generally fall into two categories: those that perform no feature\nselection (like RCTD), treating all genes equally. This often leads to the dilution of informative signals\nby noisy or ubiquitously expressed genes, reducing overall deconvolution accuracy, especially in\ncomplex or heterogeneous tissues. Another group of methods, including Spotiphy and SPOTlight,\nrelies on a top-feature selection strategy by prioritizing the most spatially informative genes. While\nthis approach enriches for highly specific markers, it inevitably discards a substantial amount of\nauxiliary transcriptomic information that may also contribute to accurate cell-type decomposition.\nIn contrast, UCASpatial employs a weighted aggregation framework that incorporates a broad set\nof genes while adaptively emphasizing those with stronger cell-type specificity. Specifically, using\nShannon entropy-based gene weighting allows UCASpatial to quantify and highlight features that\nshow informative expression patterns across cell subsets without excluding genes outright. This\nenables the model to retain a richer and more complementary transcriptomic representation, rather than\n20\nrestricting inference to a limited feature subset. Consequently, UCASpatial robustly captures subtle\ncell-type differences, reduces bias introduced by gene filtering, and maintains stable performance\nacross diverse simulated and real-world datasets.\nTo further demonstrate that the performance improvement of UCASpatial stems from its\nintegrated strategy of weighting and feature extraction optimization, we refer the reviewer to the\nsystematic ablation study we conducted in our original submission. This study involved constructing\na variant model, UCASpatial_ablation, which shares the same overall framework but lacks the core\ncomponents of entropy-based weighting and meta-purity filtering. The benchmarking clearly revealed\nthat removing these key optimizations caused the model’s performance to drop from the top rank to a\nlevel comparable with lower-performing methods (Figure 9, Figure 4, and Figure 5 for Reviewer #3).\nThis result clearly verifies that the entropy-based weighting and meta-purity filtering process are\nfundamental to UCASpatial’s superior performance over other algorithms. To ensure this crucial\nempirical evidence receives the necessary prominence, we have formally integrated these results into\nthe revised manuscript as fig. S1E-G, I, K-M, O-U, and Y-Z, fig. S2F-H and L-N, fig. S3I-L, fig.\nS4E-H, and fig. S6G-I.\nWe have expanded the discussion in the revised manuscript on Page 15, line 637, to clarify these\nmethodological distinctions and included corresponding ablation study results and additional\nbenchmarks. We sincerely thank the reviewer for encouraging us to elaborate on these aspects, which\ngreatly strengthen the interpretation and contextualization of our findings.\n(A) Graphic model showing the classification of 750 simulated spots into three groups, each representing a different\ncomplexity level - low, medium, and high, respectively. The classification is based on the entropy of cell-type\ncomposition for each spot, with each group containing 250 spots.\n(B-D) Boxplot showing spot-level root mean square error (RMSE) (B), PCC (C), and F1 score (D) for each method\nacross three complexity levels. Red boxes label the performance of UCASpatial_ablation. Boxplot elements are\ndefined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points,\noutliers. Statistical analysis was performed using a two-tailed paired t-test. ns, non-significant, p > 0.05; *0.01 ≤ p <\n0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.\n(E) Line charts showing cell-subpopulation-level Pearson correlation coefficient (PCC) for each method across 20\ncell subpopulations. ‘avg’ and ‘maxm’ represent the average and maximum PCC, respectively, of the seven methods\nexcluding UCASpatial. The red arrow labels the performance of UCASpatial_ablation.\n(F) Left, line charts representing cell-subpopulation-level PCC of CD8+ effector memory T cells (CD8+ Tem), CD8+\nexhausted T cells (CD8+ Tex), CD8+ tissue-resident memory T cells (CD8+ Trm), and natural killer cells (NK) for\n22\neach method. Right, stacked bar plot showing the accumulated count of false predictions of these cell subpopulations\nfor each method. The red arrows label the performance of UCASpatial_ablation. Colors represent different cell\nsubpopulations. FP, false positive; FN, false negative; FP + FN, total error frequency.\n(G) The cell-subpopulation-level PCC values for each method were binned into four confidence groups: 0-0.25 (low),\n0.25-0.5 (medium-low), 0.5-0.75 (medium-high), and 0.75-1 (high). Stacked bar plot showing the ratio of four\nconfidence groups for each method. The red arrow labels the performance of UCASpatial_ablation.\n(H) Graphic model representing the classification based on the proportion range of a specific cell subpopulation\nwithin the simulated spots. All simulated spots were classified into ten groups from proportion 0-10% to proportion\n0-100% for each cell subpopulation.\n(I) Line charts representing the average RMSE (left), PCC (middle), and F1 score (right) of 20 cell subpopulations\nacross various proportions for each method. Each dot in the line charts represents the average value from multiple\nsimulations. The red arrows label the performance of UCASpatial_ablation.\n3. Line 791 and Table S1, Simulation strategy: Despite adding mIHC validation, the benchmarking\nstill relies heavily on simulated data. However, the simulation strategy may not capture the full\ncomplexity of real spatial transcriptomics data, such as the spatial correlation structures between\nneighboring spots and cell-type-specific capture efficiency variations. The author should consider\nthese features to construct more reliable synthetic spots.\nResponse: We appreciate the reviewer's concern that synthetic data may fail to fully capture the\ncomplexity of real spatial transcriptomics data, particularly regarding spatial correlation structures and\ncell-type-specific capture efficiency variations.\nTo address this and construct a more reliable benchmark, we adopted a data-driven strategy\nderived from high-definition spatial transcriptomics. We generated pseudo-spots from three 10x\nVisium HD CRC slides5 that intrinsically retain local neighborhood correlations and spot-level capture\nvariability (Figure 10A for Reviewer #3). To mitigate the ambiguity that an 8-µm bin may contain\neither multiple cells or fractions of a single cell, we re-ran Space Ranger with cell segmentation and\nannotated each segmented cell (Figures 10B-C for Reviewer #3). We then partitioned each tissue\nsection into square bins of 49 µm × 49 µm—closely matching the size of a standard Visium spot—and\naggregated the transcriptomes of all segmented cells within each bin to produce a pseudo-spots\nexpression profile (Figure 10A for Reviewer #3). Critically, every pseudo-spot retains its original\nspatial coordinates and is accompanied by ground-truth cell-type proportions, providing a definitive\nreference for validation. In total, we generated 39,813 pseudo-spots across three colorectal cancer\ntissue sections, collectively preserving the intricate spatial heterogeneity and architectural complexity\nof the tumor microenvironment.\nBenchmarking on this dataset yielded results consistent with our previous findings: UCASpatial\nachieved the highest Pearson correlation coefficient (mean PCC = 0.87) and F1 score (mean F1 score\n= 0.82), as well as the lowest root mean square error (mean RMSE = 0.05) among contemporary\nmethods at the spot level (Figure 10D for Reviewer #3), even under scenarios of varying complexity\n(Figure 10E for Reviewer #3). In addition, we evaluated the performance of each method in\ndeciphering different cell subpopulations, and found that UCASpatial also performed best across PCC,\nRMSE, and F1 score (Figure 10F and G for Reviewer #3). For example, in the P1CRC section, only\nUCASpatial and stereoscope accurately estimated the proportion of CD8 Tem in CRC (Figure 10H for\nReviewer #3). However, stereoscope incorrectly annotated the proportion of B cells, whereas\nwere unable to apply cell2location to this specific analysis due to its known substantial computational\nburden.\nAdditionally, results for P2CRC and P5CRC are available in the Zenodo repository:\nhttps://zenodo.org/records/12634581?preview=1&token=eyJhbGciOiJIUzUxMiJ9.eyJpZCI6ImU0OGM4YjBj\nLTliYzQtNDc4ZC05Njg3LTczYTdhOTExNTliNiIsImRhdGEiOnt9LCJyYW5kb20iOiI0ZjYyYzU1OGYwNz\nc4M2E2NDQ5YWFiMjZkZmY2OGE4NCJ9.qcl1aa2yo6mmdkZXdHw2exPmfhyA7M5oPS1ZoKXbqVpBd\nQbSWkagYIGL2U3Dhj2ffSde0rhe5MFiuIRYuPC_Vw, due to space constraints.\nWe believe this validation using Visium HD-derived ground truth fully addresses the reviewer's\nconcern regarding the realism of our benchmarking data. We have updated the benchmarking section\nin the revised manuscript to include these results on Page 6, lines 233-246, and fig. S2I-N.\n24\n25\n26\n(A) Graphic model showing the generation of stimulated spatial transcriptomic data with ground truth.\n(B) Spatial mapping of CRC samples using single cell–resolution bins, colored by cell subpopulation annotation.\n(C) Bubble heatmap showing the expression of marker genes of 34 cell subpopulations in three Visium HD sections.\nDot size indicates the fraction of cells per subpopulation that express a given marker gene. Colors show the average\nexpression of each gene in each subpopulation.\n(D) Boxplots showing the performances PCC, RMSE, and F1 score of UCASpatial and other methods, including\nUCASpatial without entropy weighting and meta-purity filtering (UCASpatial_ablation), at spot-level. Colors\nrepresent different methods.\n(E) Boxplot showing PCC, RMSE, and F1 score for each method across three complexity levels.\n(F) Boxplot showing the capacity for distinguishing diverse cell subpopulations of UCASpatial,\nUCASpatial_ablation, and other methods. PCC, RMSE, and F1 score between the predicted proportions and the\nground truth were calculated for each cell subpopulation.\n(G) Line charts showing cell subpopulation-level PCC (left), RMSE (middle), and F1 score (right) for each method\nacross 34 cell subpopulations. ‘avg’ represents the average RMSE, PCC or F1 score, and ‘best’ represents the\nminimum RMSE, or the maximum PCC or F1 score, respectively, of the five methods excluding UCASpatial and\nUCASpatial_ablation.\n28\n(H) The proportion of cell subpopulations in the simulated spatial transcriptomics data (based on the P1CRC section)\nevaluated by each deconvolution method, with ground-truth proportion shown alongside for reference.\nBoxplot elements are defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5×\ninterquartile range; points, outliers. Statistical analysis was performed using a two-tailed paired t-test. ns, non-\nsignificant, p > 0.05; *0.01 ≤ p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.\nRemarks on code availability:\nThe UCASpatial GitHub repository provides a README file with installation instructions and sample\ncode examples for both Python and R implementations. I successfully tested the newly added Python\npackage (developed after the first review round) using the Spatiotemporal Cardiac Dynamics Post-\nInfarction Dataset (GSE163129 and GSE165857). The authors have provided reproduction code for\nFigures 2-5 in the main manuscript. However, a notable limitation for reproducibility is the lack of\nsource data for the simulation dataset and h5ad or Seurat objects for any real datasets used in the\nmanuscript. They are not available in either the data availability section of the paper or the GitHub\nrepository, preventing reproduction of any results. To improve reproducibility and community\naccessibility, I strongly recommend that the authors: (1) upload all Seurat or h5ad data files used in\nFigures 2-5, (2) provide a representative sample dataset for users to test the software, and (3) develop\na simple tutorial to facilitate easier implementation.\nResponse: We thank the reviewer for carefully evaluating the reproducibility of UCASpatial and for\nthe constructive suggestions. We fully agree that providing all source data and usage examples is\nessential for transparent and reproducible research. In response to the reviewer’s concerns, we have\nsubstantially improved the data availability and documentation resources.\nFirst, we have uploaded the code and associated Seurat data corresponding to Figures 2-6 to\nZenodo for convenient access and reproducibility. These data enable full reproduction of every analysis\npresented in the main figures of the manuscript. The datasets can be found here:\nhttps://zenodo.org/records/12634581?preview=1&token=eyJhbGciOiJIUzUxMiJ9.eyJpZCI6ImU0OGM4YjBj\nLTliYzQtNDc4ZC05Njg3LTczYTdhOTExNTliNiIsImRhdGEiOnt9LCJyYW5kb20iOiI0ZjYyYzU1OGYwNz\nc4M2E2NDQ5YWFiMjZkZmY2OGE4NCJ9.qcl1aa2yo6mmdkZXdHw2exPmfhyA7M5oPS1ZoKXbqVpBd\nQbSWkagYIGL2U3Dhj2ffSde0rhe5MFiuIRYuPC_Vw.\nSecond, to facilitate users in testing the software without requiring large dataset downloads, we\nhave added a representative example dataset and developed a comprehensive tutorial within the\nGitHub repository (see\nhttps://github.com/BIGHanLab/UCASpatial/blob/main/docs/documentation/04_UCASpatial_Examp\nle.md) as well as our dedicated tutorial website (https://bighanlab.github.io/UCASpatial/). The tutorial\nprovides detailed, step-by-step instructions using the example data (also included in the Zenodo upload\nas Example_Data_Human_CRC_ST_data.rds and Example_Data_sc_ref_downsampled.rds),\nenabling users to easily implement and test UCASpatial.\nThird, we have constructed a step-by-step tutorial website\n(https://bighanlab.github.io/UCASpatial/) (Figure 11 for Reviewer #3) that integrates installation\nguidance, usage examples for both Python and R implementations, application programming interfaces\n(API) documentation, and frequently asked questions (FAQ), further supporting community adoption.\nMoreover, we have updated the Data and Code Availability section in the revised manuscript to\nclearly state: “The data and code to reproduce the results in this study are available in the Zenodo\nrepositories.” We believe these updates significantly enhance reproducibility and accessibility, and we\n29\nthank the reviewer again for the helpful suggestions that motivated these improvements.\nFigure 11 for Reviewer #3. The tutorial website for UCASpatial (https://bighanlab.github.io/UCASpatial/).\nWe would like to express our sincere gratitude to Reviewer #3 for their insightful critique, which\nhas been instrumental in elevating the quality and robustness of our work. Your constructive feedback\ncatalyzed significant improvements across the manuscript, including the empirical validation of our\nentropy-based weighting strategy, the expansion of benchmarking to include Slide-seq V2\nhippocampus and high-resolution Visium HD data, and the development of the new Reference\nDiagnosis Tool and a dedicated tutorial website. We believe these revisions, driven by your valuable\nsuggestions, have not only strengthened the methodological rigor of our study but also substantially\nenhanced the practical utility and generalizability of UCASpatial for the broader research community.\nReferences\n1. Buechler, M.B., Pradhan, R.N., Krishnamurty, A.T., Cox, C., Calviello, A.K., Wang, A.W., Yang, Y.A., Tam, L.,\nCaothien, R., Roose-Girma, M., et al. (2021). Cross-tissue organization of the fibroblast lineage. Nature 593, 575–\n579. 10.1038/s41586-021-03549-5.\n2. Franzen, L., Olsson Lindvall, M., Huhn, M., Ptasinski, V., Setyo, L., Keith, B.P., Collin, A., Oag, S., Volckaert, T.,\nBorde, A., et al. (2024). Mapping spatially resolved transcriptomes in human and mouse pulmonary fibrosis. Nat\nGenet 56, 1725–1736. 10.1038/s41588-024-01819-2.\n3. Abedini, A., Levinsohn, J., Klotzer, K.A., Dumoulin, B., Ma, Z., Frederick, J., Dhillon, P., Balzer, M.S., Shrestha,\nR., Liu, H., et al. (2024). Single-cell multi-omic and spatial profiling of human kidneys implicates the fibrotic\nmicroenvironment in kidney disease progression. Nat Genet 56, 1712–1724. 10.1038/s41588-024-01802-x.\n4. Ma, Y., and Zhou, X. (2022). Spatially informed cell-type deconvolution for spatial transcriptomics. Nat\nBiotechnol 40, 1349–1359. 10.1038/s41587-022-01273-7.\n5. Oliveira, M.F., Romero, J.P., Chung, M., Williams, S.R., Gottscho, A.D., Gupta, A., Pilipauskas, S.E., Mohabbat,\nS., Raman, N., Sukovich, D.J., et al. (2025). High-definition spatial transcriptomic profiling of immune cell\npopulations in colorectal cancer. Nat Genet 57, 1512–1523. 10.1038/s41588-025-02193-3.\n30"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for addressing my comments.", "reviewer_3": "The manuscript has been significantly strengthened by the empirical validation of the entropy-based weighting strategy and the expanded benchmarking, which now includes Slide-seq V2 hippocampus and high-resolution Visium HD data. Furthermore, the development of the Reference Diagnosis Tool and the accompanying tutorial website effectively addresses previous concerns regarding methodology and reproducibility.\nI only have one minor comment: In the rebuttal, the authors mention the development of a \"Reference Diagnosis Tool\" (illustrated in Figure 2 for reviewer #3 of the response) to address concerns regarding reference quality and dissociation biases. While I have verified the tool's presence on the package website, I was unable to locate a specific description of this tool or the corresponding figure within the revised manuscript text or the Supplementary Materials. Please ensure that this section is explicitly included in the final manuscript or supplementary documents, as it is a valuable addition for users.\nOverall, the manuscript is significantly improved. I have no further concerns."}, "author_response": "Point-by-point response\nReviewers' Comments:\nReviewer #1:\nI thank the authors for addressing my comments.\nResponse: We sincerely thank Reviewer #1 for their time, valuable comments, and positive assessment\nof our revised manuscript. We are pleased that the reviewer is satisfied with our responses and the\nsubsequent revisions made to the work.\nReviewer #3:\nThe manuscript has been significantly strengthened by the empirical validation of the entropy-based\nweighting strategy and the expanded benchmarking, which now includes Slide-seq V2 hippocampus\nand high-resolution Visium HD data. Furthermore, the development of the Reference Diagnosis Tool\nand the accompanying tutorial website effectively addresses previous concerns regarding methodology\nand reproducibility.\nI only have one minor comment: In the rebuttal, the authors mention the development of a \"Reference\nDiagnosis Tool\" (illustrated in Figure 2 for reviewer #3 of the response) to address concerns regarding\nreference quality and dissociation biases. While I have verified the tool's presence on the package\nwebsite, I was unable to locate a specific description of this tool or the corresponding figure within\nthe revised manuscript text or the Supplementary Materials. Please ensure that this section is explicitly\nincluded in the final manuscript or supplementary documents, as it is a valuable addition for users.\nOverall, the manuscript is significantly improved. I have no further concerns.\nRemarks on code availability:\nThe development of the Reference Diagnosis Tool and the accompanying tutorial website addresses\nprevious concerns regarding methodology and reproducibility.\nResponse: We would like to express our sincere gratitude to Reviewer #3 for the positive evaluation\nand for recognizing the improvements made to our manuscript. We are particularly encouraged that\nthe new benchmarking on Slide-seq V2 and Visium HD data, as well as the reference diagnosis tool,\nmet the reviewer's expectations.\nRegarding the minor comment, we appreciate the reviewer highlighting this point. We have now\nensured that the reference diagnosis tool is formally integrated into the manuscript on the\nImplementation of the reference diagnosis tool in the Methods section.\n1"}]

Decision letter