HRCHY-CytoCommunity identifies hierarchical tissue organization in cell-type spatial maps

Domain: chemistry
Published: 2026-02-28
Editorial outcome: major_revision

Abstract
Tissues are organized through the assembly of diverse cell types into multicellular structures that exhibit hierarchical spatial organization. We present HRCHY-CytoCommunity, a graph neural network framework for identifying multi-level tissue structures directly from cell-type annotated spatial maps. It integrates differentiable graph pooling, adaptive edge pruning, and consistency and balance regularization in an end-to-end model, simultaneously inferring robust structures across multiple scales while preserving complete cellular coverage and fully nested relationships. The framework also supports cross-sample hierarchy alignment via cell-type enrichment-based clustering. Benchmarking on diverse spatial omics datasets, HRCHY-CytoCommunity outperforms existing hierarchical and non-hierarchical methods in identifying both coarse-grained tissue compartments and fine-grained cellular neighborhoods. Applied to a breast cancer cohort with clinical outcomes, the framework enables hierarchical prognostic stratification of patients and reveals survival-associated spatial patterns. HRCHY-CytoCommunity represents a general and scalable tool for deciphering tissue organization from single cells to multicellular modules, and ultimately to intact tissues and organs.

Introduction
The rapid advancement of spatial omics technologies has significantly advanced our understanding of tissue spatial organization 1. High-throughput, sequencing-based spatial transcriptomics technologies, such as Visium 2,3, Stereo-seq4, and Slide-seq 5,6, provide transcriptome-wide gene expression profiling while preserving spatial context. Imaging-based spatial transcriptomics and proteomics technologies like MERFISH7, CODEX 8, IMC9, and MIBI-TOF10, enable highly multiplexed, single-cell-resolution measurement of targeted RNA or protein markers. To systematically decipher functional tissue organization based on these spatial omics data, researchers have introduced the concepts of cellular neighborhoods and spatial domains, which serve as multicellular modules where distinct cell types interact and collaborate to support tissue functions 8,11,12. This has spurred the development of numerous computational methods for identifying such structural modules 13–24.

Many tissues exhibit inherent hierarchical organization, wherein large structures contain smaller, functionally specialized substructures. This pattern is evident in the layered architecture of the brain, the zonated morphology of lymphoid organs such as the spleen, and the compartmentalized ecosystem of tumors 25–27. For example, some tumor tissues can be coarsely partitioned into neoplastic, immune, and stromal compartments 28,29, which can be further subdivided into finer neighborhoods enriched with specific cell types, such as B cells and T cells, granulocytes, or fibroblasts 30. These hierarchical patterns carry biological and clinical significance. As an example for illustrating the role of coarse-grained tissue structures, in triple-negative breast cancer (TNBC), spatial segregation between neoplastic and immune compartments is associated with favorable outcomes 10. As a representative fine-grained multicellular module, the presence of B and T cell-enriched tertiary lymphoid structure serves as a positive prognostic marker across multiple cancer types 31–33. Uncovering such multi-level organization is therefore essential for under standing how tissues assemble from individual cells into functional modules.

Despite its biological importance, few computational methods are designed specifically to identify hierarchical tissue structures. The current state-of-the-art method, NeST 34, identifies nested hierarchies in spatial transcriptomics data by detecting co-expression hotspots, which are defined as groups of spatially colocalized cells that co-express a set of genes. However, NeST relies on gene expression features, limiting its performance on data with sparse molecular measurements, such as spatial proteomics datasets 8,9. Moreover, the structures identified by NeST often fail to cover all cells in the sample and may lack clearly nested relationship across scales. For instance, cells within the same fine-grained structure can be assigned to different coarse-grained structures, complicating biological interpretation. Additionally, NeST treats spatially disconnected hotspots as distinct tissue structures, potentially impeding the discovery of structures with spatially discontinuous distribution.

Non-hierarchical methods 15,20–23, though not originally designed for hierarchical analysis, can be repurposed to identify multi-scale tissue structures through post hoc adjustment of clustering resolution. However, this strategy requires each hierarchical level to be trained and identified independently, often resulting in a lack of nested relationships across scales. Moreover, the absence of joint optimization or constraints between coarse- and fine-grained representations may compromise biological interpretability and reduce the reliability of the resulting hierarchical structures.

To overcome these limitations, we developed HRCHY-CytoCommunity, a hierarchical extension of the CytoCommunity framework 16, that identifies multi-scale tissue organization in cell-type spatial maps. Unlike conventional multi-step approaches, HRCHY‑CytoCommunity employs a unified end-to-end learning framework that simultaneously infers tissue structures at multiple levels while preserving hierarchical consistency through joint optimization. By leveraging differentiable graph pooling, adaptive edge pruning, and consistency and balance regularization during training, the model identifies robust and fully nested tissue structures that encompass all cells and reflect biologically meaningful spatial hierarchies. This joint learning strategy enables bidirectional information flow across scales, promoting the identification of structurally consistent and functionally relevant tissue architectures. The model also supports cross-sample integration through cell-type enrichment-based clustering, aligning hierarchical structures across samples for comparative analysis.

In this study, we focused on a two-level hierarchy, defining coarse-grained structures as tissue compartments (TCs) and fine-grained structures as cellular neighborhoods (CNs). We demonstrated the performance of HRCHY-CytoCommunity across a wide range of simulated and real spatial omics datasets spanning diverse tissues, technologies, and modalities. Benchmarking results showed that our method outperformed existing hierarchical and non-hierarchical approaches in reconstructing biologically meaningful spatial hierarchies. Finally, using a clinical breast cancer cohort, we illustrated how these hierarchical structures can serve as multi-scale prognostic indicators, highlighting the potential of HRCHY-CytoCommunity to support translational research.

Results
Overview of HRCHY-CytoCommunity
HRCHY-CytoCommunity is an unsupervised graph neural network framework designed to identify hierarchical multicellular structures from single-cell spatial maps that integrate both cell type and spatial location information. The framework comprises a base module for soft hierarchical tissue structure assignment and a consistency and balance regularization module to ensure robust identification. Fig. 1 illustrates the identification process of tissue structures at two hierarchical levels.

In the base module (Fig. 1a; Methods), HRCHY-CytoCommunity begins by constructing a K-nearest-neighbor (KNN)-based cell-cell proximity graph, where nodes represent cells and edges connect spatially adjacent cells. Each node is equipped with an attribute vector encoding cell type information. The model then applies a graph convolution layer with ReLU activation to generate node embeddings, which are subsequently transformed through a fully-connected layer with Softmax activation into soft assignment vectors for fine-grained tissue structures (i.e., CNs). These vectors represent the probability of each cell belonging to a specific CN.

HRCHY-CytoCommunity applies a differentiable graph pooling layer to coarsen the original graph into a new completed graph, where each pooled node corresponds to a fine-grained CN. The edge weights in this coarsened graph reflect inter-CN connectivity strength. To prevent over-smoothing effect from message passing during graph convolution, an adaptive edge-pruning step is applied on the coarsened graph before performing a second graph convolution and fully-connected layer to produce soft assignment vectors for coarse-grained tissue structures (i.e., TCs). These vectors represent the probability of each CN belonging to a specific TC. The learning process is guided by two graph minimum cut (MinCut)-based loss functions 35, separately optimizing fine-grained CN and coarse-grained TC assignments.

The second module (Fig. 1b; Methods) introduces both consistency and balance regularization to enhance the accuracy and robustness of hierarchical tissue structure identification. During training, multiple perturbed versions of the cell‑cell proximity graph are generated by randomly dropping node features (DropNode). These perturbed graphs as fed into the base module to produce soft fine-grained CN and coarse-grained TC assignment matrices. Consistency regularization is applied to enforce agreement among the resulting assignments from different perturbed graphs, while balance regularization based on entropy prevents cluster collapse. The total loss function combines the base MinCut loss with the two regularization terms. Final robust hierarchical structures are obtained by performing hard assignment on the stable soft assignment matrices.

Since HRCHY-CytoCommunity operates as an unsupervised learning framework optimized for single-sample analysis, an optional cross-sample tissue structure alignment module is included (Supplementary Fig. 1; Methods). This module computes cell-type enrichment scores for each identified tissue structure across all samples, concatenates them into a composite feature matrix, and applies clustering to align structures into a unified set. This allows the identification of shared architectural patterns while preserving sample‑specific multicellular structures.

For comprehensive evaluation, we compared HRCHY-CytoCommunity not only with NeST 34, a method specifically designed for identifying hierarchical tissue structures, but also with six widely-used non-hierarchical tissue structure identification methods. This comparison was conducted to assess whether explicit hierarchical modeling yields more biologically meaningful results than repeatedly adjusting the resolution of flat structural identification methods to infer multi-level organization. The methods include five deep learning-based approaches (CellCharter 21, NicheCompass 20, GraphST 15, SpaGCN 22, and SpaSEG 23) as well as the seminal spatial domain detection method (Giotto Suite 24) (Supplementary Data 1).

Performance evaluation on coarse-grained TC identification
To evaluate the ability of HRCHY-CytoCommunity in identifying coarse-grained tissue structures, we first applied it to an imaging-based spatial proteomics dataset of mouse spleen generated through the Co-Detection by Indexing (CODEX) technology 31. This dataset comprises three healthy mouse spleen samples (BALBc-1, BALBc-2, and BALBc-3) with an average of 81,760 cells per sample. Each sample was annotated into 27 distinct cell types (Fig. 2a, left column) and two major functional compartments: the red pulp and the lymphoid compartment, the latter comprising the periarteriolar lymphoid sheath, B-zone, and marginal zone 36–39 . These manually annotated compartments served as the gold standard for the coarse-grained TCs (Fig. 2a, right column). To ensure a rigorous and fair comparison, HRCHY-CytoCommunity was evaluated against NeST and four non-hierarchical tissue structure identification methods (Supplementary Data 1). However, Giotto Suite was excluded from this analysis due to its computational limitations in handling large-scale spatial omics datasets. SpaSEG was also excluded from benchmarking in this dataset because it is not applicable for spatial proteomics data (Supplementary Note 7).

Using a cluster stability-based criterion to determine the optimal numbers of coarse-grained TCs and fine-grained CNs (Methods and Supplementary Fig. 2a), HRCHY-CytoCommunity identified two TCs and seven CNs. Across all samples, the two TCs corresponded accurately to the red pulp and the lymphoid compartment. Notably, although the lymphoid compartment is a spatially discontinuous structure, HRCHY-CytoCommunity successfully captured its distribution (Fig. 2b, left column). In contrast, NeST identified nested hierarchical structures starting from single-gene hotspots 34, thus some cells were not assigned to any TC (labeled as unidentified). Moreover, NeST often split the spatially discontinuous lymphoid compartment into multiple distinct TCs (Fig. 2b, right column). Non-hierarchical tissue structure identification methods also exhibited limited accuracy in delineating coarse-grained TCs. For example, CellCharter captured only B cell-enriched regions within the lymphoid compartment (Fig. 2c, first column), while NicheCompass and SpaGCN detected only T cell (CD4+ and CD8+)-enriched areas of the lymphoid compartment (Fig. 2c, second and fourth columns). GraphST correctly identified TC distributions in sample BALBc-1 but failed to distinguish the lymphoid compartment in the other two samples (Fig. 2c, third column). To quantitatively evaluate performance, we measured the concordance between identified TCs and manually annotated compartments using adjusted mutual information (AMI) and Macro-F1 scores (Methods). The results demonstrated that HRCHY-CytoCommunity achieved significantly higher AMI and Macro-F1 scores than most state-of-the-art methods (one-sided paired t-test P-values < 0.05), while exhibiting performance comparable to GraphST (Fig. 2d).

To further evaluate the ability of HRCHY-CytoCommunity to identify coarse-grained tissue structures in sequencing-based spatial omics data, we applied it to a spatial transcriptomic dataset of human colorectal cancer (CRC) generated through the Visium HD technology 3, which profiles 18,085 genes (Supplementary Data 1). This dataset comprises three samples (named as P1CRC, P2CRC, and P5CRC). On average, each sample contained 420,502 bins at 8 μm resolution covering 10 major cell types (Fig. 2e, left column) and two major compartments: the tumor and the normal tissue compartment, serving as the gold standard for coarse-grained TCs (Fig. 2e, right column). For benchmarking, HRCHY-CytoCommunity was compared with NeST and three non-hierarchical spatial domain detection methods (Supplementary Data 1). Due to memory constraints associated with processing large-scale spatial omics data, GraphST, SpaGCN and Giotto Suite were excluded from the analysis (Supplementary Note 7).

Based on the cluster stability criterion, HRCHY-CytoCommunity identified two coarse-grained TCs and ten fine-grained CNs (Supplementary Fig. 3a). In sample P2CRC, HRCHY-CytoCommunity accurately distinguished tumor and normal tissue compartments. Similarly, in samples P1CRC and P5CRC, the two compartments were clearly delineated, though with moderate inclusion of intestinal epithelial-enriched regions within the tumor compartment (Fig. 2f, left column). In contrast, NeST misclassified the majority of spatial bins as tumor regions and left most normal tissue regions unassigned (labeled as unidentified), limiting the interpretability of the resulting tissue structures (Fig. 2f, right column). Among the non-hierarchical methods, CellCharter exhibited competitive overall performance relative to HRCHY-CytoCommunity, but it misclassified epithelial-enriched regions as tumor compartments in sample P2CRC (Fig. 2g, first column). Meanwhile, NicheCompass incorrectly assigned large numbers of immune and stromal cells to the tumor compartment, while SpaSEG accurately identified the tumor compartments but split the normal compartments into multiple distinct TCs (Fig. 2g, second and third columns). Quantitative evaluation confirmed that HRCHY-CytoCommunity outperformed other methods, showing a statistically significant improvement in Macro-F1 score (one-sided paired t-test, P-values < 0.05) and a performance level comparable to CellCharter. On the AMI metric, HRCHY-CytoCommunity achieved a higher average score, albeit without statistical significance (Fig. 2h).

We further investigated the biological relevance of the fine-grained CNs identified by HRCHY-CytoCommunity in these two datasets. In the absence of ground-truth labels for fine-grained structures, HRCHY-CytoCommunity consistently uncovered spatially coherent and biologically interpretable CNs that aligned with established anatomical and functional structures (Supplementary Note 1). For instance, in the mouse spleen, CNs accurately corresponded to anatomical regions such as the B-cell zone and the periarteriolar lymphoid sheath, while in the human CRC, CNs represented distinct tumor-core and tumor-edge regions along with immune cell infiltrations.

In summary, HRCHY-CytoCommunity outperforms existing state-of-the-art hierarchical and non-hierarchical tissue structure identification methods in identifying coarse-grained TCs across evaluated datasets. Moreover, the fine-grained CNs detected by our method show stronger correspondence to known biological structures, demonstrating its ability to reveal meaningful multi-scale spatial organizations in complex tissues.

Performance evaluation on fine-grained CN identification
To further assess the capability of HRCHY-CytoCommunity in identifying fine-grained CNs, we first applied it to an imaging-based spatial transcriptomics dataset of the healthy mouse hypothalamic preoptic region. The dataset was generated using the Multiplexed Error-Robust Fluorescence in situ Hybridization (MERFISH) technology 7, profiling 155 genes. It comprises five samples representing distinct brain regions, named as Bregma-0.04, Bregma-0.14, Bregma+0.06, Bregma+0.16, and Bregma+0.26 based on their relative distances to the bregma. Each sample contains multiple small and symmetric hypothalamic nuclei regions, which were delineated in the original study through manual inspection of well-characterized brain histology 7,40. We previously manually mapped these nuclei onto single-cell spatial maps that served as the ground truth for fine-grained CN evaluation (Fig.3a). For benchmarking, HRCHY-CytoCommunity was compared with NeST and six non-hierarchical spatial domain detection methods (Supplementary Data 1).

Using the cluster stability criterion, HRCHY-CytoCommunity identified 12 fine-grained CNs and two coarse-grained TCs (Supplementary Fig. 4a). This model successfully identified multiple symmetric CNs that corresponded to established hypothalamic nuclei (Fig.3b, left column). For instance, BNST (Bed Nucleus of the Stria Terminalis), LPO (Lateral Preoptic Area) and MPA (Medial Preoptic Area) were consistently identified across all five samples. Relatively larger nuclei, such as ACA (Anterior Commissure, Anterior Part) and MPN (Medial Preoptic Nucleus), were detected in most samples, while relatively smaller nuclei, including PS (Parastrial Nucleus), Pe (Periventricular Nucleus), and PaAP (Paraventricular Nucleus) were identified in several samples. Our method also captured sample-specific structures, such as PVA (Paraventricular Nucleus, Anterior Part) and Fx (Fornix) in the Bregma-0.14 sample and SHy (Suprachiasmatic nucleus shell) in the Bregma+0.06 sample. In contrast, NeST identified only one or two CNs per sample (Fig.3b, right column), showing poor correspondence to the manual annotations. This limitation likely stems from NeST’s reliance on gene expression features without pre-specifying the number of tissue structures and inability to identify spatially discontinuous parts of the same nuclear structure. Comparison with six additional non-hierarchical spatial domain detection methods (Supplementary Data 1) revealed that while most methods preserved general symmetry of identified CNs, their results often contained intermixed CNs, leading to suboptimal performance (Fig. 3c). Among them, GraphST and CellCharter performed relatively well. GraphST correctly identified relatively large CNs such as BNST, but intermixed CN boundaries reduced its accuracy. CellCharter achieved superior identification of BNST, MPA, and MPN in one sample (Bregma-0.14). However, these methods consistently failed to capture smaller regions like PS and PaAP across other samples, resulting in overall inferior performance compared to HRCHY-CytoCommunity. Quantitative evaluation demonstrated that HRCHY-CytoCommunity significantly outperformed NeST (one-sided paired t-test, P-values < 0.001) and all other non-hierarchical methods (one-sided paired t-test, P-values < 0.05) in both AMI and Macro-F1 metrics, with the exception of CellCharter on the AMI score and GraphST on the Macro-F1 score, where the differences were not statistically significant (Fig. 3d).

To further evaluate the performance of HRCHY-CytoCommunity on sequencing-based spatial transcriptomics data, we applied it to a mouse intracerebral hemorrhage Stereo-seq dataset, which profiled 23,625 genes. This analysis focused on nine representative samples, with an average of 81,831 cells per sample, encompassing 27 molecularly defined regions that served as the ground truth for CNs (Supplementary Fig. 5a). HRCHY-CytoCommunity was compared with NeST and four non-hierarchical spatial domain detection methods for benchmarking (Supplementary Data 1). Giotto Suite and SpaSEG were excluded due to excessive memory requirements (Supplementary Note 7).

HRCHY-CytoCommunity identified 20 fine-grained CNs and three coarse-grained TCs (Supplementary Fig. 6a), successfully reconstructing symmetric CNs throughout the mouse brain. Specifically, our method accurately resolved biological meaningful CNs including cortical layers 1 to 6, olfactory areas (OLF), meninges, and pallidum. In contrast, CNs identified by NeST cannot correspond to known molecularly defined regions in most samples (Supplementary Fig. 5b). Among the non-hierarchical methods, while NicheCompass produced the most competitive results benefiting from its cell-cell communication modeling, HRCHY-CytoCommunity demonstrated superior sensitivity in detecting subtle structures such as OLF in the Naive-1 sample and islands of Calleja in the D1-3 sample (Supplementary Fig. 5b, left column), which were not well captured by NicheCompass (Supplementary Fig. 5c, second column). This enhanced sensitivity may arise from HRCHY-CytoCommunity’s use of cell-type annotations as input features, enabling better discrimination of closely related cell subtypes, though this design occasionally resulted in over-segmentation. Quantitatively, HRCHY-CytoCommunity achieved significantly higher AMI and Macro-F1 scores than NeST (one-sided paired t-test, P-values < 1E-7) and other non-hierarchical methods (one-sided paired t-test, P-values < 0.001), with performance comparable to NicheCompass (Supplementary Fig. 5d).

We further examined the coarse-grained TCs identified by HRCHY-CytoCommunity in these two datasets. The analysis revealed that our method consistently detected biologically meaningful TCs across multiple samples, including distinct compartments enriched for neuronal subtypes (e.g., excitatory and inhibitory neurons) and non-neuronal cells in the hypothalamic data, as well as anatomical regions such as the brain stem, cortex and striatum in the intracerebral hemorrhage data (Supplementary Note 2; Supplementary Figs. 4c and 6c). These findings confirm that HRCHY-CytoCommunity robustly reconstructs hierarchically organized tissue structures that correspond to biologically relevant domains, demonstrating its applicability across both imaging-based and sequencing-based spatial transcriptomics platforms.

Cross-platform generalization of HRCHY-CytoCommunity to low-resolution spatial transcriptomics
To evaluate the generalization ability of HRCHY-CytoCommunity on spot-based spatial transcriptomics data lacking single-cell-resolution cell-type annotations, we applied it to two representative datasets: a mouse hippocampus Slide-seq V2 dataset and a human breast cancer Visium V1 dataset (Supplementary Data 1). For each dataset, we first inferred the cell-type composition of each spot using RCTD 41, a state-of-the-art deconvolution method (Supplementary Figs. 7a and 8a). The inferred cell-type fractions were used as node attributes to construct a spot-spot proximity graph, which was then input into HRCHY-CytoCommunity.

We first evaluated the method on the mouse hippocampus Slide-seq V2 data (10 μm-spot resolution), using the annotated Allen Mouse Brain Atlas 42 as the ground-truth reference (Figs. 4a and b, left panels). HRCHY-CytoCommunity identified three coarse-grained TCs and 11 fine-grained CNs (Supplementary Fig. 7b), showing strong alignment with the reference atlas. Specifically, at the coarse-grained level, the identified TCs accurately corresponded to major anatomical structures, including the cortical region, hippocampus, and brain stem (Fig. 4a, middle panel). At the fine-grained level, the method successfully recovered substructures, such as CA1-CA3 regions, dentate gyrus, corpus callosum, and distinct cortical layers (Fig. 4b, middle panel). In contrast, although NeST detected several CNs that matched the reference atlas (e.g., CA1-CA3 regions, dentate gyrus, and corpus callosum), it produced small and fragmented substructures rather than large-scale coherent TCs at the coarse-grained level. Moreover, NeST left a considerable proportion of spots unassigned (labeled as unidentified) for both TCs and CNs, limiting the interpretability of the hierarchical tissue organization (Figs. 4a and b, right panels).

We next applied HRCHY-CytoCommunity to the human breast cancer Visium V1 data (55 μm-spot resolution), using manual annotations of hierarchical tissue structures from Xu et al.43 as ground-truth reference (Figs. 4c and d, left panels). HRCHY-CytoCommunity identified two coarse-grained TCs and nine fine-grained CNs (Supplementary Fig. 8b). The method accurately distinguished tumor and non-tumor compartments (Fig. 4c, middle panel), whereas NeST incorrectly partitioned discontinuous tumor regions into multiple fragmented TCs (Fig. 4c, right panel). At the fine-grained level, both methods identified CNs consistent with manual annotations, such as DCIS/LCIS_1, DCIS/LCIS_3, and IDC_4, but NeST left the Healthy_1 region unassigned, reducing its interpretability (Fig. 4d, middle and right panels). Notably, both methods reclassified the manually annotated IDNC_3 tumor region as a CN (e.g., CN-8 in our results) belonging to the non-tumor compartment (Fig. 4d, red arrowhead). Cell-type enrichment analysis revealed significant co-enrichment of B-cells and T-cells within this CN (Methods; one-sided Wilcoxon rank-sum test, P-values < 0.001; Fig. 4e). Differential expression analysis (|log₂FC| > 0.25) further indicated strong humoral immune activity in this region, characterized by up-regulation of immunoglobulin constant region genes (e.g., IGHG1/3/4, IGHA1, IGHM, IGKC, IGLC1/2/3 and JCHATIN), consistent with a plasma cell-enriched microenvironment 44–46 (Fig. 4f). These results suggest that the original manual annotation of IDC_3 as tumor tissue was likely inaccurate, and that HRCHY-CytoCommunity more reliably captured the underlying biologically relevant tissue organization.

In summary, HRCHY-CytoCommunity generalizes effectively across low-resolution spatial transcriptomics platforms, consistently identifying biologically meaningful hierarchical tissue structures even in the absence of high-quality single-cell-resolution cell-type annotations. In contrast to NeST, which often produces fragmented substructures and leaves many spots unassigned, our method captures both large- and small-scale spatially coherent tissue structures and achieves complete cellular coverage.

Cross-sample hierarchical tissue structures in triple-negative breast cancer
To further evaluate whether the hierarchical tissue structures identified by HRCHY-CytoCommunity can capture inter-patient heterogeneity in multicellular organization within tumor tissues, we applied our method to a triple-negative breast cancer (TNBC) dataset generated using the Multiplexed Ion Beam Imaging by Time-Of-Flight (MIBI-TOF) technology 10 (Supplementary Data 1). This dataset consists of 15 patients with compartmentalized tumors, characterized by clearly spatial segregation between immune and neoplastic cell populations. Our analysis had two main objectives: (1) to assess whether HRCHY-CytoCommunity-identified coarse-grained TCs accurately correspond to immune and neoplastic cell-dominated regions compared to NeST; and (2) to investigate whether fine-grained CNs reflect inter-patient heterogeneity even among tumors broadly classified as compartmentalized at the macroscopic level.

HRCHY-CytoCommunity identified two TCs and 14 fine-grained CNs across the dataset (Supplementary Fig. 9a). In most samples, our method successfully delineated one TC dominated by immune cells and another by neoplastic cells, consistent with the spatial distribution patterns of cell types expected in compartmentalized tumors. In contrast, NeST identified only a single TC in multiple samples (e.g., patients 9, 10 and 16), likely due to limitations in feature availability preventing simultaneous identification of both compartments (Fig. 5a; Supplementary Fig. 9b). In patient 28, where immune and neoplastic cells are spatially contiguous, both methods effectively identified the TCs. However, in patient 32, where the tumor exhibits spatially discontinuous immune and neoplastic cell populations, HRCHY-CytoCommunity correctly separated them into distinct TCs, whereas NeST partitioned them into five distinct TCs (Fig. 5a, fourth row). To quantitatively evaluate TC identification performance, we calculated the proportion of neoplastic and immune cells that were correctly assigned to neoplastic and immune cell-dominated TCs. HRCHY-CytoCommunity showed significantly better performance across the 15 patients than NeST (one-sided paired t-test, P-value = 0.024; Fig. 5b).

While compartmentalized tumors can be broadly categorized into immune and neoplastic TCs, fine-grained CN analysis can reveal deeper heterogeneity within these tumors. We performed cross-sample integration analysis of CNs from all patients and generated a unified set of 26 CNs (Methods; Fig.5c; Supplementary Fig. 10). We found that these unified CNs captured both conserved and patient-specific biological features. For example, both patients 4 and 9 contained unified CN-6 (neoplastic cell-enriched) in their neoplastic TCs and unified CN-22 (CD4+ T cell- and B cell-enriched) in their immune TCs. However, unified CN-8 (macrophage-, CD8+ T cell- and endothelial cell-enriched) was unique to patient 4’s immune TC, and unified CN-20 (neoplastic-neutrophil mixed) appeared only in patient 9’s neoplastic TC (Figs. 5c and d). These findings indicate that even among compartmentalized tumors, individuals maintain distinct multicellular organizational patterns.

We further stratified the patients by performing hierarchical clustering based on the unified CN composition profiles of each sample. While some patients (e.g., patients 3 and 32) exhibited highly similar unified CN compositions, others (e.g., patients 9 and 10) showed clear divergence (Fig.5e, left; Supplementary Fig. 10). To determine whether these clustering patterns reflect inherent spatial organizational heterogeneity rather than mere differences in cell-type abundance, we also performed the clustering using cell-type composition profiles. The results indicated that although some patients with similar cell-type compositions (e.g., patients 28 and 35) also shared comparable CN compositions, it was more common to observe marked differences in unified CN configurations among patients with highly similar cell-type distributions, as exemplified by patients 4 and 41 (Fig.5e; Supplementary Fig. 10). These findings suggest that fine-grained CNs capture inter-patient heterogeneity in spatial tissue organization that is not fully explained by cell-type proportion alone. Therefore, while coarse-grained TC analysis can distinguish compartmentalized from mixed tumors 10, the analysis of fine-grained CNs offers potential for further subclassification of patients.

Hierarchical tissue structure-based stratification of breast cancer patients
To validate the biological significance and clinical relevance of the hierarchical tissue structures identified by HRCHY-CytoCommunity, we applied the method to a breast cancer spatial proteomics dataset generated using the imaging mass cytometry (IMC) technology 9. The dataset consists of 263 breast cancer patients, including 190 survivors and 73 deceased individuals. HRCHY-CytoCommunity identified 30 unified TCs and 50 unified CNs across the cohort.

We first stratified patients into five distinct groups by performing hierarchical clustering based on the presence or absence of unified TCs as features (Fig. 6a, top) and conducted survival analysis using associated clinical data (Fig. 6b). Patients in Group 3, characterized by the presence of TC-13, exhibited significantly better prognosis compared to the remaining cohorts (hazard ratio = 0.2, P-value = 0.031; Fig. 6c). TNC-13 represents a coarse-grained structure significantly enriched with CK7+ neoplastic cells (Fig. 6a, bottom), consistent with previous findings that CK7+ neoplastic cell phenotype is associated with favorable clinical outcomes 7.

We next focused on Group 1, which contained the largest number of patients but lacked obvious shared coarse‑grained TC characteristics. To investigate whether finer‑grained CN features could further stratify this group, we used the presence or absence of a specific unified CN as the clustering feature. This analysis revealed CNs with significant prognostic value (Figs. 6d and e; Supplementary Fig. 11), and we explored their potential functions through cell‑type enrichment analysis (Methods; Fig. 6f). For example, within Group 1, a subgroup of 29 patients containing unified CN-21 showed significantly worse prognosis (hazard ratio = 2.2, P-value = 0.026; Fig. 6d). CN-21 was significantly enriched with small elongated fibroblasts in all 29 patients, and in some cases also exhibited significant enrichment of small circular fibroblasts or diverse cancer-associated fibroblast (CAF) subtypes, including Fibronectinhi, Vimentinhi, and SMAhi Vimentinhi fibroblasts (Fig. 6f). Accumulating evidence indicates that CAFs promote tumor growth, metastatic and therapeutic resistance 47–49, and may form a dense cellular barrier that impedes immune-neoplastic cell interactions 48,50. In this study, we observed that CN‑21, enriched with normal fibroblasts and CAFs, functioned as a fence‑like structure separating immune cell‑enriched CNs (e.g., CN‑6 and CN‑40) from neoplastic cell‑enriched CNs (e.g., CN‑8, CN‑36, and CN‑28), as illustrated in patients P14, P23, and P112 (Fig. 6g; Supplementary Fig. 12).

Another notable example was a 13-patient subgroup within Group 1 containing unified CN-34, which was associated with a markedly elevated risk of death (hazard ratio = 5.4, P-value < 0.001; Fig. 6e). CN-34 was significantly enriched with T cell subpopulation-1 in all 13 patients, and in some patients also showed significant enrichment of macrophages (Fig. 6f). This CN exhibited a spatially dispersed distribution surrounded by diverse neoplastic cell‑ or CAF‑enriched CNs (e.g., CN‑1, CN‑3, and CN‑47), as observed in patients P54, P80, and P191 (Fig. 6h; Supplementary Fig. 12). This spatial pattern suggests that immune cell infiltration into the tumor core may be obstructed, potentially leading to regional immune evasion 51.

In summary, these findings demonstrate that HRCHY-CytoCommunity enables hierarchical patient stratification by leveraging multi-scale tissue structures from coarse-grained TCs to fine-grained CNs. This approach reveals patient subgroups with distinct survival implications across different spatial scales of tissue organization.

Ablation, robustness, and scalability analyses of HRCHY-CytoCommunity
To systematically evaluate the stability and scalability of HRCHY-CytoCommunity, we performed a series of complementary analyses, including ablation studies, sensitivity analyses, robustness assessments and scalability analyses (Supplementary Notes 3-6).

Comprehensive ablation experiments demonstrated that each key component of the HRCHY-CytoCommunity model, including consistency and balance regularization, adaptive edge pruning, and adaptive α scheduling (Methods) , contributed positively to overall performance (Supplementary Note 3; Supplementary Figs. 13-15). Sensitivity analyses indicated that the model maintained stable performance across a wide range of hyperparameter settings, including DropNode rate, λ_consis, λ_balance, the number of neighbors K in the KNN graph (cell-cell proximity graph) construction, and the number of perturbations used for generating perturbed cell-cell proximity graphs (Supplementary Note 4; Supplementary Figs. 16-18). Robustness evaluation against cell-type label inaccuracies and variations in annotation resolution further confirmed that the model consistently identified reliable hierarchical tissue structures under noisy or heterogeneous input annotations (Supplementary Note 5; Supplementary Figs. 19 and 20). Finally, scalability analyses on large-scale datasets (containing up to 473k spots) demonstrated high computational efficiency in both runtime and memory usage. For example, on a dataset of 473k spots, HRCHY-CytoCommunity completed hierarchical tissue structure identification within 8.02 minutes, using only 8.69 GB of memory, confirming its capability to handle large-scale spatial omics datasets (Supplementary Note 6; Supplementary Figs. 21 and 22).

Taken together, these results demonstrate that HRCHY-CytoCommunity is robust, scalable, and generalizable across diverse data modalities.

Methods
Soft hierarchical tissue structure assignment
We began by constructing an undirected KNN graph as the cell-cell proximity graph to represent the spatial omics data, where each node corresponds to a cell or spot. For single-cell-resolution datasets, the node attribute vector was constructed using one-hot encoding to capture cell type information (Fig. 1a). For low-resolution datasets, the node attribute vector was constructed using the estimated cell-type composition of each spot. The graph was built by computing Euclidean distances between nodes (cells) based on their spatial coordinates, connecting each node to its K nearest neighbors (excluding itself).

For the undirected KNN graph with n nodes, we applied a single graph convolution layer 54 with the ReLU activation function to produce a node embedding matrix 𝑍(1) ∈ ℝ𝑛×𝑑 ,
formulated as
𝑍(1) = ReLU(GNN1(𝑋(1),𝐴(1);𝜃GNN1)) (1)
where each row of 𝑍(1) is a learned d-dimensional embedding vector of a cell. 𝑋(1) ∈ ℝ𝑛×𝑡 denotes the initial cell-type attribute matrix, and t is the total number of cell types. 𝐴(1) ∈ {0,1}𝑛×𝑛 is the adjacency matrix of the undirected KNN graph. The graph convolution operator was defined as
𝐳(𝟏) = ReLU(𝜃1 𝐱(𝟏) +𝜃2 ∑ 𝐱(𝟏) ) (2)
𝐢 𝐣
𝑗∈𝑁(𝑖)
where 𝐱 is the embedding vector of node 𝑖, and 𝑁(𝑖) denotes the first-order neighborhood derived from the matrix 𝐴(1). 𝜃 and 𝜃 are trainable parameters in the graph neural network GNN . The embedding dimension d was empirically set to 128 in this study.
1
Next, we used a fully-connected neural network with no hidden layers (also referred to as a linear layer) and a Softmax activation function to convert the node embedding matrix 𝑍(1) ∈ ℝ𝑛×𝑑 into a soft fine-grained tissue structure assignment matrix 𝑆(1) ∈ ℝ𝑛×𝑐1, which was formulated as below.
𝑆(1) = Softmax(FC1(𝑍(1);𝜃FC1)) (3)
where each element of 𝑆(1) represents the probability of a cell (row) belonging to one of the c1 fine-grained tissue structures (column). 𝜃 represents trainable parameters in the fully-connected neural network FC . The hyperparameter c specifies the maximum number of fine-grained tissue structures to be detected and the model automatically determines the final number (less than or equal to c1).
1
To further identify coarse-grained tissue structures, we applied a differentiable graph pooling layer 35,55 to generate a coarsened graph with c1 pooled nodes, where each pooled node represents a fine-grained tissue structure. The adjacency matrix 𝐴(2) ∈ ℝ𝑐1×𝑐1 and node embedding matrix 𝑋(2) ∈ ℝ𝑐1×𝑑 of the coarsened graph were formulated as follows.
𝐴(2) = 𝑆(1)𝑇 𝐴(1)𝑆(1) (4)
𝑋(2) = 𝑆(1) 𝑇 𝑍(1) (5)
The resulting coarsened graph is fully-connected with self-loops. Self-loops were removed and the adjacency matrix was normalized as
𝐴̂ = 𝐴(2) −𝐼 𝑑𝑖𝑎𝑔(𝐴(2)) (6)
𝑐1
𝐴𝑝𝑜𝑜𝑙 = 𝐷̂− 2𝐴̂𝐷̂− 2 (7)
where 𝐷̂ is the degree matrix of 𝐴̂. Each entry in 𝐴𝑝𝑜𝑜𝑙 represents the normalized strength of the connection between two fine-grained tissue structures.

To relieve the smoothing effect of message-passing operation, we introduced an adaptive edge-pruning threshold 𝑇 . Only edges with weights in the adjacency matrix 𝐴𝑝𝑜𝑜𝑙 greater than this threshold are retained and their weights are reset to 1.
𝐴̂(2) = { 1 𝑖𝑓 𝐴𝑝𝑜𝑜𝑙 𝑖𝑗 > 𝑇 (8)
0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Since excessive or insufficient edge pruning may affect model performance (Supplementary Note 3; Supplementary Figs. 13b-15b), the adaptive threshold was defined as
𝑇 = 1 (9)
𝑐1−1
This value corresponds to the case where all fine-grained tissue structures are connected with equal strength.

Finally, another graph convolution and fully-connected layers were applied to generate the soft coarse-grained tissue structure assignment matrix:
𝑍(2) = ReLU(GNN2(𝑋(2),𝐴̂(2);𝜃GNN2)) (10)
𝑆(2) = Softmax(FC2(𝑍(2);𝜃FC2)) (11)
Each element of 𝑆(2) ∈ ℝ𝑐1×𝑐2 represents the probability of a fine-grained structure (row) belonging to one of the c2 coarse-grained structures (column). 𝑍(2) ∈ ℝ𝑐1×𝑑 is the updated embedding matrix of pooled nodes in the coarsened graph. 𝜃 and 𝜃 represent trainable parameters in the graph neural network GNN and fully-connected neural network FC , respectively. The hyperparameter c specifies the maximum number of coarse-grained tissue structures to be detected and the model automatically determines the final number (less than or equal to c2).
2
The loss function with regard to the soft hierarchical tissue structure assignment (base module) was defined as follows.
𝐿 = 𝛼×𝐿 +(1−𝛼)×𝐿 (12)
𝐵𝑎𝑠𝑒 𝐹𝑖𝑛𝑒 𝐶𝑜𝑎𝑟𝑠𝑒
where 𝛼 is a weight parameter used for balancing the fine-grained loss 𝐿 and the coarse-grained loss 𝐿 . We further introduced an adaptive scheduling strategy for determining 𝛼 to dynamically balance these two losses throughout training. Specifically, 𝛼 was scheduled to decay from 0.9 to 0.1 during training, shifting focus from fine-grained to coarse-grained structure learning. Both 𝐿 and 𝐿 use the graph MinCut-based formulation, optimizing the matrix 𝑆(1) and 𝑆(2), respectively.
𝐿 = − ∑ 𝑗 𝑐 = 1 1 (𝑆(1)𝑇 𝐴(1)𝑆(1)) 𝑗𝑗 +‖ 𝑆(1)𝑇 𝑆(1) − 𝐼𝑐1‖ (13)
𝐹𝑖𝑛𝑒 ∑ 𝑐
𝑗=1
1
(𝑆(1)𝑇 𝐷(1)𝑆(1)) 𝑗𝑗
‖𝑆(1)𝑇 𝑆(1)‖
𝐹
√𝑐1
𝐹
𝐿 = − ∑ 𝑐 𝑘 = 1 2 (𝑆(2) 𝑇 𝐴(2)𝑆(2)) 𝑘𝑘 +‖ 𝑆(2) 𝑇 𝑆(2) − 𝐼𝑐2‖ (14)
𝐶𝑜𝑎𝑟𝑠𝑒 ∑ 𝑐
𝑘=1
2
(𝑆(2)𝑇 𝐷(2)𝑆(2)) 𝑘𝑘
‖𝑆(2)𝑇 𝑆(2)‖
𝐹
√𝑐2
𝐹
where 𝐷(1) ∈ ℝ𝑛×𝑛 and 𝐷(2) ∈ ℝ𝑐1×𝑐1 are degree matrices derived from the undirected KNN graph and the coarsened graph, respectively. Both 𝐿 and 𝐿 consist of two terms. The first term encourages clustering of strongly connected nodes, and the second term encourages orthogonal and balanced tissue structure assignments 30.

Robust hierarchical tissue structure identification
To enhance the stability of HRCHY-CytoCommunity, we introduced a consistency and balance regularization module for identifying robust hierarchical tissue structures (Fig. 1b). Specifically, we adapted the GRAND framework 56 to the hierarchical tissue structure assignment task. The core idea involves stochastically generating multiple perturbed versions of the cell-cell proximity graph (i.e., KNN graph) by randomly dropping node features, and then enforcing cluster assignment consistency across these perturbed graphs. An entropy-based balance regularization was also used to prevent cluster collapse. This module improves both the robustness and accuracy of hierarchical tissue structure identification (Supplementary Note 3; Supplementary Figs. 13a-15a).

First, given the original KNN graph with adjacency matrix 𝐴(1) and node attribute matrix 𝑋(1), we randomly removed the entire attribute vectors of a subset of nodes. This dropout strategy was referred to as DropNode. In contrast to standard dropout, which only masks individual elements of 𝑋(1), DropNode removes all features of selected nodes, thereby accounting for graph structural effects and explicitly reducing reliance on specific neighboring nodes. This generates more stochastic perturbations and has been shown to achieve better performance than standard dropout 56.

Formally, for each node 𝑖, a binary mask 𝜖 was sampled from a Bernoulli distribution 𝜖 ~Bernoulli(1−𝛿), where 𝛿 is a predefined DropNode rate. The perturbed feature vector was then computed as 𝐱̃ = 𝜖 ∙𝐱 , where 𝐱 represents the original feature vector of node 𝑖. The perturbed feature matrix 𝑋̃ was used only during training. During inference, features were rescaled by multiplying (1−𝛿) to match the expected activation magnitude of 𝑋̃.

Next, the soft hierarchical tissue structure assignment module was performed on each perturbed cell-cell proximity graph. By repeating the DropNode procedure M times, M sets of predicted labels for both fine-grained CNs and coarse-grained TCs were obtained. To evaluate prediction consistency, we calculated the normalized assignment distribution center for each node:
𝐬̅ = 1 ∑𝑀 𝐬 (𝐦) (15)
𝐢 𝑚 𝑚=1 𝐢
where 𝐬 denotes the predicted tissue structure assignment distribution of node 𝑖 in the 𝑚-th perturbed graph.

Finally, the consistency regularization loss was defined as the average 𝐿 distance between the assignment distributions from each perturbed graph and the distribution center:
𝐿 = 1 ∑𝑀 ∑𝑛 ‖𝐬 (𝐦) −𝐬̅ ‖ (16)
𝐶𝑜𝑛𝑠𝑖𝑠 𝑀 𝑚=1 𝑖=1 𝐢 𝐢 2
This loss encourages prediction consistency under multiple stochastic perturbations, thereby improving the robustness of hierarchical tissue structure assignments produced by HRCHY-CytoCommunity.

However, consistency regularization alone may increase the risk of the collapse of the coarse-grained assignment layer (i.e., all nodes being assigned to a single cluster), trivially minimizing the consistency loss. To prevent such collapse, we introduced an entropy-based balance regularization term:
𝐿 = 1− 1 ∑𝑀 𝐇 (2) (17)
𝐵𝑎𝑙𝑎𝑛𝑐𝑒 𝑚=1
𝑀 log𝑐2
𝐇 (2) = −∑𝑐2 𝑠 (2) log𝑠 (2) (18)
𝐦 𝑘=1 𝑚,𝑘 𝑚,𝑘
where 𝐇 is the entropy of the cluster assignment distribution for the m-th perturbation, and c2 is the number of coarse-grained TCs. This term encourages balanced cluster assignments by maximizing the entropy of the predicted distributions across clusters.

The overall loss function of HRCHY-CytoCommunity was then defined as
𝐿 = 𝐿 +𝜆 ×𝐿 +𝜆 ×𝐿 (19)
𝑇𝑜𝑡𝑎𝑙 𝐵𝑎𝑠𝑒 1 𝐶𝑜𝑛𝑠𝑖𝑠 2 𝐵𝑎𝑙𝑎𝑛𝑐𝑒
where 𝐿 is the base loss from the soft hierarchical tissue structure assignment, 𝐿 is the consistency regularization loss, and 𝐿 is the entropy-based balance regularization loss. The hyperparameters 𝜆 and 𝜆 control the relative contributions of these terms. A comprehensive list of the hyperparameter settings used for all datasets was provided in Supplementary Table 1.

Determination of optimal numbers of hierarchical tissue structures based on cluster stability
To identify hierarchical tissue structures using HRCHY-CytoCommunity, users are required to pre-specify the desired numbers of fine-grained CNs, denoted as c₁, and coarse-grained TCs, denoted as c₂. Inspired by the approach implemented in CellCharter, HRCHY-CytoCommunity employs a cluster stability-based procedure to determine the optimal values of the pair (c₁, c₂).
Specifically, the procedure begins by defining a search range for both c₁ and c₂, as well as the number of independent runs R. For each candidate (c₁, c₂) pair, HRCHY-CytoCommunity performed R runs. To evaluate the consistency of the resulting clusters (i.e., CNs and TCs) across runs, we computed the Fowlkes-Mallows Index (FMI) 57, which quantifies the similarity between two cluster assignments as follows:
𝐹𝑀𝐼 = 𝑇𝑃 (20)
√(𝑇𝑃+𝐹𝑃)(𝑇𝑃+𝐹𝑁)
where TP, FP, and FN represent true positives, false positives, and false negatives, respectively, evaluated over pairs of cells assigned to the same or different clusters. The FMI ranges from 0 to 1, with higher values indicating greater similarity between cluster assignments.

In HRCHY-CytoCommunity, the FMI was computed separately for the identified CNs and TCs across all pairs of runs for a given (c₁, c₂) configuration. This yielded a set of CN-level FMIs and a set of TC-level FMIs. For each level, the median FMI across all run pairs was computed to summarize the cluster stability at that level. The final stability score for a (c₁, c₂) configuration was then defined as the average of the CN-level and TC-level median FMIs. The optimal (c₁, c₂) pair was selected as the one that yields the highest stability score. Among the R models associated with the optimal (c₁, c₂) configuration, the model achieving the lowest training loss was chosen to output the final hierarchical tissue structures.

In this study, for datasets where the ground-truth number of CNs or TCs is known, we fixed the known value and applied the cluster stability-based procedure to determine the optimal number of structures at the other hierarchical level. For multi-sample datasets, the optimal (c₁, c₂) pair was selected based on the highest average cluster stability score across all samples.

Unified hierarchical tissue structure generation across samples
To facilitate integrated analysis of multiple tissue samples, we developed a clustering-based module to align hierarchical tissue structures across different samples (Supplementary Fig. 1).
First, cell-type enrichment scores were computed for each tissue structure (including both coarse-grained and fine-grained structures) within every sample using the hypergeometric test. Then, we horizontally concatenated the cell-type enrichment score matrices from all samples to construct a composite feature matrix. This composite matrix was then input into a clustering algorithm to group tissue structures from different samples into a unified set of hierarchical tissue structures.

For the human TNBC MIBI-TOF dataset and the human breast cancer IMC dataset, we applied hierarchical clustering with Ward’s method. Euclidean distance was used to measure the similarity between patients’ tissue structure profiles.

Running of published methods
We compared the performance of HRCHY-CytoCommunity with seven established tissue structure identification methods, including six widely-used non-hierarchical methods, GraphST 15, NicheCompass 20, CellCharter 21, SpaGCN 22, Giotto Suite 24, and SpaSEG 23, as well as NeST 34, a state-of-the-art method specifically designed for identifying nested hierarchical structures in spatial transcriptomics data. To ensure a fair and consistent comparison, the number of TCs and CNs to be detected by each method was set to match the ground-truth annotations provided in the original studies.

The Python code base for NeST was downloaded from https://github.com/bwalker1/NeST. This code base was applied to seven datasets (Supplementary Data 1). In accordance with NeST requirements, we utilized protein or mRNA expression data and cell spatial coordinates as inputs. For benchmarking purposes, we considered the first layer of the NeST output as coarse-grained TCs and the last layer as the fine-grained CNs. In cases where a cell population was assigned to multiple tissue structures within the same layer simultaneously, cells were allocated to the tissue structure with the purpose of obtaining highest accuracy. For imaging-based datasets with limited numbers of features, the model hyperparameters were adjusted to ensure the identification of multiple single-gene hotspots and hierarchical structures. Notably, NeST does not require pre-specification of the number of tissue structures but demands the configuration of eight other hyperparameters. Detailed settings were listed in Supplementary Table 2.

The Python package GraphST (v.1.1.1) was applied to three datasets (Supplementary Data 1). For all three datasets, default hyperparameters were used for model training, following the official tutorial at https://deepst-tutorials.readthedocs.io/en/latest/. Notably, for the mouse spleen CODEX dataset, all 30 protein markers were used as highly variable genes (HVGs) due to the limited features. Spatial clustering was subsequently performed using the mclust package in R.

The Python package NicheCompass (v.0.3.0) was applied to four datasets (Supplementary Data 1). For the mouse spleen CODEX dataset, we followed the single-sample tutorial at https://nichecompass.readthedocs.io/en/latest/tutorials/notebooks/mouse_cns_single_sample.html for better performance. For the other three datasets, we followed the sample-integration tutorial at https://nichecompass.readthedocs.io/en/latest/tutorials/notebooks/mouse_cns_sample_integration.html. Due to the large numbers of cells in the mouse ICH Stereo-seq dataset and the human CRC Visium HD dataset, edge_batch_size was set to 1024. For the mouse spleen CODEX dataset, the scanpy.pp.neighbors and scanpy.tl.louvain functions were used with resolution = 0.005, 0.002, and 0.005 to identify coarse-grained TCs and 0.055, 0.035, and 0.045 to identify fine-grained CNs for samples BALBc-1, BALBc-2, and BALBc-3, respectively. For the mouse hypothalamic preoptic region MERFISH dataset, the parameter resolution was set to 0.02 and 0.6 for all five samples to identify coarse-grained TCs and fine-grained CNs, respectively. For the mouse ICH Stereo-seq dataset, the parameter resolution was set to 0.05 and 0.8 for all eight samples to identify coarse-grained TCs and fine-grained CNs, respectively. For the human CRC Visium HD dataset, the parameter resolution was set to 0.009 and 0.1 for all three samples to identify coarse-grained TCs and fine-grained CNs, respectively.

The Python package CellCharter (v.0.3.5) was applied to four datasets (Supplementary Data 1). Following the official tutorial at https://cellcharter.readthedocs.io/en/latest/index.html, scArches 58 was used for dimension reduction and integration for the mouse spleen CODEX dataset. For other spatial transcriptomics datasets, scVI 59 was used for dimension reduction and integration. Default hyperparameters were used to retrieve aggregated cell embeddings. The gmm.predict function was used to identify TCs and CNs.

The Python package SpaGCN (v.1.2.7) was applied to three datasets (Supplementary Data 1). All parameters were set according to the official tutorial at https://github.com/jianhuupenn/SpaGCN/blob/master/tutorial/tutorial.ipynb, except that max_epochs was increased to 200 for training convergence. Notably, for the mouse spleen CODEX dataset, the number of principal components was set to 30 due to limited numbers of features.

The Python implementation of SpaSEG was obtained from https://github.com/y-bai/SpaSEG. SpaSEG was applied to two datasets (Supplementary Data 1). All parameters were set following the official tutorial at https://github.com/y-bai/SpaSEG/blob/main/DLPFC_151673_demo.ipynb. During data preprocessing, the appropriate platforms were selected according to the datasets’ platforms.

The R package Giotto Suite (v4.2.1) was applied to the MERFISH dataset following the two-stage workflow recommended in the official tutorial. First, genes with significant spatial patterns were identified using the spatial co-expression module detection method at https://giottosuite.com/articles/spatial_coexpression_modules.html. Then the identified genes were used as features for spatial clustering with the HMRF (Hidden Markov Random Field) model 60(https://giottosuite.com/articles/hmrf.html).

Cell-type deconvolution using RCTD
Cell-type composition of each spot in the mouse hippocampus Slide-seq V2 dataset and the human breast cancer Visium V1 dataset was inferred using the RCTD 41 algorithm implemented in the R package spacexr (v2.2.0). Following the official tutorial (https://github.com/dmcable/spacexr ), we set CELL_MIN_INSTANCE = 1 and UMI_min = 1 in the create.RCTD function for both datasets. In the run.RCTD function, we set doublet_mode = multi for the Visium V1 dataset, and set doublet_mode = doublet for the Slide-seq V2 dataset, based on the different spatial resolutions of these two technologies. For each dataset, single-cell RNA-seq data from matched tissue types were used as reference to train the RCTD model (Data Availability). The resulting cell-type proportion matrices were used as the input node attributes in HRCHY-CytoCommunity.

Quantitative performance evaluation using CODEX, MERFISH, Stereo-seq and Visium HD datasets
For the mouse spleen CODEX dataset, the ground-truth assignments of cells to fine-grained CNs, including red pulp, marginal zone, B-cell zone and periarteriolar lymphoid sheath (PALS), were derived from the original study 31. For the mouse hypothalamic preoptic region MERFISH dataset, the outlines of hypothalamic nuclei regions were obtained from the original study 5, and the ground-truth assignments of cells to CNs were manually annotated based on these outlines. For the mouse ICH Stereo-seq datatset, the CN annotations were obtained from the original study61. For the human CRC Visium HD dataset, the TC annotations were obtained from the original study 3. We quantitatively evaluated the performance of HRCHY-CytoCommunity and benchmarked methods using two evaluation metrics: Macro-F1 score and Adjusted Mutual Information (AMI) with calculations performed using the Python package scikit-learn (v1.1.2). The formulas for these evaluation metrics are as follows.
𝐹1 𝑠𝑐𝑜𝑟𝑒 = 2×(𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛×𝑅𝑒𝑐𝑎𝑙𝑙) (21)
𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛+𝑅𝑒𝑐𝑎𝑙𝑙
𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛 = 𝑇𝑃 (22)
𝑇𝑃+𝐹𝑃
𝑅𝑒𝑐𝑎𝑙𝑙 = 𝑇𝑃 (23)
𝑇𝑃+𝐹𝑁
𝐴𝑀𝐼 = 𝐼(𝐺𝑇;𝑇𝑆)−𝐸{𝐼(𝐺𝑇;𝑇𝑆)} (24)
1 [𝐻(𝐺𝑇)+𝐻(𝑇𝑆)]−𝐸{𝐼(𝐺𝑇;𝑇𝑆)}
2
where true positives (TP), true negatives (TN), false positives (FP) and false negatives (FN) represent the comparisons between the predicted tissue structures and the ground-truth cell assignments. 𝐸{𝐼(𝐺𝑇;𝑇𝑆)} denotes the expected mutual information between the ground-truth assignments of cells (GT) and the tissue structure predictions (TS). 𝐻(𝐺𝑇) and 𝐻(𝑇𝑆) represent the entropies of the ground-truth assignments of cells and the tissue structure predictions, respectively. The Macro-F1 score was calculated as the average of the F1 scores across all types of ground-truth tissue structures within the dataset. The AMI is a metric based on Shannon information theory 62 that assesses the agreement between the predicted tissue structures and the ground-truth assignments, adjusting for the possibility of random agreement.

Cell-type enrichment scores of tissue structures
To perform a quantitative assessment of cell-type composition within tissue structures, we defined an enrichment score for each cell type in each tissue structure as -log10(P-value). For a single-cell spatial omics sample, the P-value was determined through a hypergeometric test, taking into account four key parameters: (1) the number of cells of the given type within the tissue structure, (2) the total number of cells within the tissue structure, (3) the number of cells of that type in the sample, and (4) the total number of cells in the sample.

For a low-resolution (spot-based) spatial transcriptomics sample (e.g., Visium V1), where each spot contains a mixture of cell types, enrichment analysis was performed using cell-type deconvolution results. Specifically, for a cell type in an identified tissue structure, we compared the cell-type fractions across all spots within the structure against those from all other spots in the sample using a one-sided Wilcoxon rank-sum test to determine whether this cell type was significantly enriched within the identified structure. To account for multiple comparisons, the resulting P-values were adjusted with the Benjamini-Hochberg method 63.

Survival analysis
To validate clinical relevance of the hierarchical tissue structures identified by HRCHY-CytoCommunity, we performed survival analysis using unified TCs or unified CNs as features in combination with clinical survival data from the breast cancer IMC dataset 7. All Kaplan-Meier survival curves were generated using the R package survival (v3.6-4). To evaluate the independent prognostic value of TC- or CN-derived features, we employed multivariable Cox proportional hazard regression models, adjusting for clinical covariates including patient age, tumor grade and tumor size.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "HRCHY-CytoCommunity presents a novel graph neural network (GNN)-based framework for identifying hierarchical\nmulticellular structures in spatial omics data. Leveraging differentiable graph pooling and graph pruning, the method robustly\nuncovers nested tissue compartments and cellular neighborhoods from cell-type spatial maps.\nMajor advances of this approach:\n1) HRCHY-CytoCommunity outperforms the state-of-the-art NeST method across multiple datasets (mouse spleen CODEX,\nhypothalamic MERFISH, TNBC MIBI-TOF). It accurately identifies spatially discontinuous structures and achieves higher\naccuracy, adjusted mutual information, and Macro-F1 scores.\n2) The hierarchical majority voting strategy enhances reproducibility while preserving nested relationships between TCs and\nCNs. Unlike NeST, it assigns all cells to structures, avoiding coverage gaps.\n3) A cell-type enrichment-based clustering module aligns structures across samples, revealing conserved and patient-\nspecific biological features (e.g., 26 unified CNs in TNBC).\n4) The method reveals functional hierarchies, such as neuron/non-neuron interactions in the hypothalamus and immune-\ntumor compartmentalization in TNBC, aligning with known biology. Applied to breast cancer, HRCHY-CytoCommunity\nenables multi-scale prognostic stratification.\nLimitations:\n1) Researchers need to pre-define/select structure counts/levels. Are there any quantitative indexes to guide the hyper-\nparameter settings?\n2) Can this model be adapted to RNA-seq based spot-structure data, such as Visium HD & StereoSeq?\nOverall, this work is a significant advance in spatial omics analysis with biological and clinical insights.", "reviewer_2": "In this study, Wang et al. introduced a spatial clustering framework HRCHY-CytoCommunity. Unlike conventional spatial\nclustering methods such as SpaGCN or graphST that deliver a single level spatial domain clustering, the proposed method\ncan present hierarchical spatial clustering, detecting nested spatial domains and structures that can not be detected with the\nconventional methods. Considering the complexity of the tissue microenvironment in various complex diseases such as\ncancer, the capability to detect hierarchical spatial clustering can help us better understand the spatial patterns underlying\nvarious biological conditions and diseases, therefore it surely will be of interest of an ever-growing community of single-cell\nspatial genomics. Although I am fully convinced the importance of the problem and appreciate the authors effort in\naddressing this issue, I found that the authors may need to enhance the benchmarking (current only one single method was\ncompared and benchmarked) and add a ablation study (to showcase the necessity of the all modules that they utilized).\nBeyond showing the effectiveness of the method and its superior performance in comparison with other method, the authors\nmay also want to show more evidence to justify why the proposed touchup major voting strategy (NOT a integral part of the\nmodel and thus can’t be optimized all together) is solid and won’t damage the performance dramatically compared to other\nintegrative alternative strategy (see major comment #2)\nMajor concerns\n(1) A potential limitation of the HRCHY-CytoCommunity method lies in its reliance on user-specified hyperparameters c1\nand c2, which define the number of fine-grained cellular neighborhoods (CNs) and coarse-grained tissue compartments\n(TCs), respectively. While the authors acknowledge that the model can learn fewer clusters than specified if overestimated,\nthere is no principled strategy provided for selecting appropriate values of c1 and c2 in practice. This imposes a significant\nconstraint on the method’s generalizability, particularly in real-world settings where prior knowledge about hierarchical\ntissue architecture is limited or absent. Without an automated or data-driven approach to determine optimal resolution, users\nmust rely on ad hoc or domain-specific heuristics, which may lead to suboptimal or biologically misleading interpretations. A\ndiscussion on how to systematically tune or infer these hyperparameter, ideally integrated to be part of the model, would\nsubstantially improve the applicability and robustness of the framework.\n(2) The hierarchical majority voting strategy used in HRCHY-CytoCommunity, while helpful for improving stability and\nenforcing a nested structure, is applied post hoc and is not integrated into the model’s training objective. This disconnection\nintroduces two key issues. First, the learned soft assignments are not explicitly optimized to yield stable or hierarchically\nconsistent outputs, which may lead to mismatches between training-time objectives and final outputs. Second, because the\nvoting procedure is non-differentiable, it prevents end-to-end training and may obscure the interpretation of how hierarchical\nstructures are learned. More critically, this touch-up procedure could obscure systematic instability in the underlying GNN\nassignments, which would not be evident without the ensemble step. As a result, the biological or clinical significance of the\nidentified structures may hinge on a potentially brittle post-processing routine.\nAs a practical alternative, the authors could consider incorporating a consistency regularization strategy, which encourages\nalignment of cluster assignments across stochastic runs during training, or a differentiable nesting regularization term, which\nsoftly penalizes misalignment between fine- and coarse-grained structures. These approaches would enable the model to\nlearn hierarchical relationships more robustly and transparently during optimization. One quick strategy is to introduce\nconsistency regularization during training. This involves running the model multiple times with stochastic perturbations (e.g.,\ndropout or graph augmentation) and penalizing discrepancies in the predicted cluster assignments across runs. By\nminimizing divergence (e.g., KL divergence) between outputs, the model is encouraged to produce stable and robust\nhierarchical assignments without relying on post hoc voting.\nAt a minimum, the authors should implement and compare one of these integrative alternatives against the current\nensemble-based (touchup major voting based) strategy to demonstrate that the model's performance and robustness are not\noverly dependent on post hoc correction. Such a comparison would strengthen the claims of generalizability and biological\ninterpretability. If the performances are indeed similar (even slightly worse), then the authors may not need to update all the\nresults with the new model structure, which is quite intensive.\n(3) The developed method integrates multiple components, including a hierarchical GNN with differentiable pooling, a post\nhoc majority-voting ensemble for robustness, and a clustering-based module for cross-sample alignment. However, the\nmanuscript does not provide an ablation study to evaluate the individual contribution of these modules to overall\nperformance. In particular, it remains unclear how much of the reported improvement over baseline methods arises from the\nGNN-based hierarchical learning itself, versus the ensemble voting or cross-sample unification steps. An ablation study\nremoving or modifying each module, such as training without voting, disabling edge pruning, or skipping the alignment step,\nwould be valuable for confirming that each component is essential and for understanding where the gains in performance\nand robustness originate.\n(4) The current benchmarking is restricted to a single method, NeST, which, while hierarchical, is designed for high-\ndimensional spatial transcriptomics and does not generalize well to low-dimensional or cell-type-based spatial data such as\nspatial proteomics. This narrow comparison limits the ability to assess the broader utility and competitiveness of HRCHY-\nCytoCommunity. In particular, several more recent and widely used frameworks—such as SpaSEG, Giotto, and GraphST—\nare highly relevant comparators. Although these methods do not explicitly model hierarchical structures, they can be run at\nmultiple resolutions or applied recursively to simulate hierarchical segmentation. Benchmarking against such methods\nwould provide a more comprehensive evaluation of whether HRCHY-CytoCommunity's explicit hierarchical modeling offers\npractical advantages in terms of spatial coherence, robustness, and biological relevance. The authors are strongly\nencouraged to include comparisons to these additional methods or, at a minimum, to discuss their relative performance.\n(5) The reliance on cell-type annotations as input features introduces a dependence on upstream preprocessing quality,\nwhich is not evaluated. The authors should assess how robust their method is to label inaccuracies or varying annotation\nresolutions, or clarify the assumptions about cell-type consistency across samples.\n(6) Given the multi-step pipeline and ensemble-based framework, the scalability of HRCHY-CytoCommunity in terms of\nruntime and memory is not discussed. The authors should provide runtime/memory profiling or clarify the computational\nrequirements to facilitate practical adoption.\nMinor comments\n(1) HRCHY-CytoCommunity models a fixed two-level hierarchy. This assumption may limit applicability to tissues with\ndeeper or non-uniform hierarchical organization. The authors should discuss whether and how the framework can be\nextended to more flexible or adaptive hierarchical depths.\n(2) The interpretability of the discovered tissue structures is not systematically evaluated beyond cell-type enrichment.\nIt would be valuable to assess how consistently the model discovers known anatomical or functional regions across samples\nand tissues, and how uninterpretable clusters are handled.\n(3) a 2-levels hierarchical structure”, ⇒ “a 2-level hierarchical structure”\n(4) Use of inconsistent capitalization for “Cellular Neighborhoods” and “Tissue Compartments” (sometimes capitalized,\nsometimes not). Maybe the authors can pick one consistent style (e.g., lowercase unless at the start of a sentence).\n(5) Line 179, “the coarse-grained red pulp compartment could be further divided into two fine-grained substructures, red pulp\nand marginal zone”. This phrasing is somewhat confusing, as it suggests that “red pulp” is both a coarse-grained and a fine-\ngrained region, and it categorizes the marginal zone, which is typically considered a transitional or fine-grained region\nbetween red and white pulp, as a substructure of red pulp. This may need clarification or correction, especially for readers\nwith histological expertise.", "reviewer_3": "The manuscript reviewed introduces HRCHY-CytoCommunity, a method for identifying hierarchical tissue structures using\nonly spatial graphs and cell-type annotations. The study presents a method to perform hierarchical structure identification\nusing graph neural networks with pruning and pooling. The study exemplified this strategy in multiple datasets, highlighting\nthe identification of discontinuous spatial patterns and presenting an improved performance in datasets with few expression\nfeatures (e.g., protein-based data).\nWhile these are promising goals, the manuscript relies heavily on the assumption that hierarchical organization in tissues\ncan be fully captured using cell-type annotations alone. This may be an oversimplification of complex tissues where cellular\nbehavior is driven by signaling gradients, cellular interactions, or dynamic gene programs. We believe this together with\nothers exposed below, limit the current version of the manuscript from being suitable for publication. In our opinion, only a\ndeeply revised version of the manuscript, with clarification of the following concerns should be considered for publication.\nMajor Comments\n1. Assumption of Hierarchy: The method assumes that tissue structures are organized hierarchically, which is not universally\ntrue. While this may apply to some biological systems (e.g., lymphoid tissues, brain), in others, spatial patterns are more\ncontinuous or transitional. The manuscript should explicitly discuss the biological validity of this assumption and its\nlimitations.\n2. Evaluation on More Diverse Datasets: The manuscript focuses primarily on datasets with well-separated cell types and\npre-annotated structures. To evaluate the robustness of HRCHY-CytoCommunity, we recommend including sequencing-\nbased spatial transcriptomics datasets—especially those from the NeST study or similar benchmarks—to compare with\nmethods like NeST, which rely on transcriptome-wide data.\n3. Limitations of Using Cell-Type Annotation Only: The exclusive reliance on cell-type labels poses a major limitation.\nHigh-quality annotations may not always be available, and annotation resolution varies by dataset. Furthermore, critical spatial\nstructures may be defined not by cell types but by cell-cell interactions, gradients, or activity states. These limitations should\nbe acknowledged, and potential solutions or future directions discussed. For instance, the author could compare the results\nof two different levels of annotation (finer vs coarse)..\n4. Comparison with Additional Methods: The current comparisons omit several state-of-the-art spatial domain identification\ntools, including SpaGCN, NicheCompass, CellCharter, Gaston-Mix among others. While these do not explicitly identify\nhierarchical domains, many can be run iteratively or parameter-tuned to extract structures at different scales. Including such\ncomparisons would provide a fairer performance benchmark.\n5. Loss Threshold Parameter (β): The parameter β is introduced as a threshold for selecting structures based on graph loss\nbut is not reflected in the mathematical formulations (Equations 7–9). The manuscript should be updated to incorporate β\ndirectly into the equations, with a formal explanation of its role and its effect on structure selection.\n6. Details of Model Parameters: Please add a table summarizing the hyperparameters used for each dataset: number of\nGNN layers, MLP layers, values for α, β, and pruning thresholds. Without this information, reproducibility and understanding\nof model performance across datasets are limited.\n7. Sensitivity Analysis: The study would benefit from conducting and reporting a sensitivity analysis for α, β, and edge\npruning on at least one image-based and one sequencing-based dataset. This is critical to demonstrate model stability and\nthe impact of key hyperparameters on performance.\n8. Choice and Estimation of Number of Structures (c₁, c₂): The model requires the number of fine- and coarse-grained\nstructures to be pre-defined. The manuscript should clarify whether these are learned or fixed, and how users should choose\nthem. If large values are set, how does the model handle oversegmentation? A discussion or heuristic strategy would be\nhelpful.\n9. Interpretability of Nested Structures: The manuscript suggests that assigning each fine-grained structure to a single\ncoarse-grained one improves interpretability. However, biological systems often contain overlapping or ambiguous\nstructures. For example, cells in a boundary region may legitimately belong to multiple domains. Instead of enforcing strict\nassignments, the model might benefit from allowing flexible memberships or highlighting ambiguity.\n10. Survival Analysis Details: The survival analysis lacks details on the inclusion of clinical covariates. Without covariate\nadjustment, survival differences may be confounded. If covariates (e.g., age, sex, tumor grade) were not included, either re-\nrun the analysis with appropriate adjustments or justify why this analysis is still valid.\n11. Validation with Simulated Data: Consider including a synthetic dataset with known hierarchical ground truth to\nbenchmark model accuracy and interpretability under controlled conditions.\n12. Multiple Sample Integration: The method currently runs HRCHY-CytoCommunity independently per sample, followed by\nclustering enrichment profiles. This strategy might not fully capture shared spatial structure across samples. Consider\ndiscussing or exploring alternatives, such as joint embedding approaches or graph alignment techniques.\nMinor Comments\n1. Please add a legend to Figure 3a to clarify the visualized elements and improve readability.\n2. The rationale for selecting the “3TC” annotation in the bregma mouse brain sample is not provided. Please justify its use\nor clarify its biological relevance.\n3. For the brain example, it would be informative to overlay the results with reference annotations from the Allen Brain Atlas\nto validate spatial domains at both coarse and fine levels.\n4. In some examples, fine-grained domains appear more as transitional gradients than discrete clusters. This challenges the\nassumption of discrete, nested hierarchies and should be discussed.\n5. While this method builds on Cytocommunity, the distinction lies primarily in the added CNN for hierarchical domain\ninference. Please clarify what conceptual or performance advancements are made beyond the original Cytocommunity\napproach.\n6. The claim that NeST is not applicable to low-feature datasets needs further explanation. Since NeST uses gene\nexpression, it may still be applicable depending on the resolution. If excluded from comparisons, the manuscript should\nclearly justify why.\n7. The use of manual annotations in unorganized organs (e.g., lung, kidney) is potentially problematic, as these organs lack\nstrong spatial hierarchy. Please explain the rationale behind these choices and whether such annotations are used as\ntraining or just evaluation.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to reviewers\nWe would like to sincerely thank the editor and the reviewers for their constructive comments\nand insightful suggestions, which have greatly helped us improve both the clarity and quality\nof our manuscript. In summary, we have addressed the following five aspects, which are\ndescribed in detail in the point-by-point responses to the reviewer comments below.\n1. Model modification and performance improvement\n• We replaced the post hoc majority-voting mechanism in the original version of\nHRCHY-CytoCommunity with the consistency and balance regularization terms that\nare directly incorporated into the model’s optimization objective. These modifications\nimproved both the accuracy and stability of hierarchical tissue structure identification\nand allowed fully end-to-end differentiable training (see Methods subsection “Robust\nhierarchical tissue structure identification” and Fig. 1b).\n• We introduced a cluster stability-based heuristic strategy to guide the determination of\noptimal numbers of hierarchical tissue structures, thereby enhancing the method’s\npracticality (see Methods subsection “Determination of optimal numbers of\nhierarchical tissue structures based on cluster stability”).\n2. Expanded evaluation and additional benchmarking\n• As suggested by Reviewers 2 and 3, we expanded the benchmarking to include six\nadditional state-of-the-art non-hierarchical tissue structure identification methods\n(SpaGCN, SpaSEG, Giotto Suite, GraphST, NicheCompass, and CellCharter), besides\nNeST. The results demonstrated that HRCHY-CytoCommunity achieved consistently\nsuperior or competitive performance in identifying hierarchical tissue structures (see\nFigs. 2-4, Supplementary Figs. 2-6, and Supplementary Notes 1 and 2).\n• We further applied HRCHY-CytoCommunity to four additional sequencing-based\nspatial transcriptomics datasets spanning diverse resolutions and technologies,\nconfirming its broad applicability across imaging- and sequencing-based platforms (see\nFigs. 2e-h, Fig. 4, Supplementary Figs. 3, 5-8 and Supplementary Notes 1 and 2).\n3. Comprehensive ablation, sensitivity, robustness, and scalability analyses\n• We added comprehensive ablation and sensitivity analyses to evaluate the contribution\nof each model component (i.e., consistency and balance regularization, adaptive edge\npruning, and adaptive α scheduling) and sensitivity of each key hyperparameter (i.e.,\nthe DropNode rate, λ_consis, λ_balance, the number of neighbors K in the KNN graph\n(cell-cell proximity graph) construction, and the number of perturbations used for\ngenerating perturbed cell-cell proximity graphs) (see Supplementary Figs. 13-18 and\nSupplementary Notes 3 and 4).\n• We conducted robustness assessments against cell-type label inaccuracies and\nannotation resolution variations, demonstrating that HRCHY-CytoCommunity remains\nstable and robust under noisy or heterogeneous input annotations (see Supplementary\nFigs. 19 and 20, and Supplementary Note 5).\n• We further performed scalability analysis, showing that HRCHY-CytoCommunity is\nhighly efficient in both runtime and memory usage. On a large-scale 473k-spot dataset,\nthe full model of HRCHY-CytoCommunity completed hierarchical tissue structures\nidentification in 8.02 minutes, requiring only 8.69 GB of memory, confirming its strong\nscalability to large-scale spatial omics datasets (see Supplementary Figs. 21 and 22,\nand Supplementary Note 6).\n4. Enhancement in survival analysis and biological interpretation\n• All survival analyses were revised to use multivariable Cox proportional hazards models\ninstead of log-rank tests, allowing adjustment for clinical covariates and providing more\nreliable effect estimates. (see Methods subsection “Survival analysis” and Figs. 6c-e).\n• We expanded applications of HRCHY-CytoCommunity across multiple tissues (e.g.,\nbrain, colorectal cancer, and breast cancer) and provided clear biological interpretations\nof the identified coarse- and fine-grained structures based on their histological contexts.\n5. Clarifications, textual revisions, and presentation improvements\n• We clarified the methodological scope of the manuscript, refined figure legends and\ndescriptions for improved readability, and ensured that all key parameters and equations\n(e.g., the reformulated loss functions) are explicitly defined.\nOverall, besides modifications and expansions of existing analyses and original figures, we\nhave added a new main figure (Figure 4), 20 new supplementary figures, nine new\nsupplementary notes, and four new supplementary tables. With these additional analyses\nmotivated by the editor’s and reviewers’ insightful comments, we believe that the manuscript\nhas been substantially strengthened.\n- Figure 4. Evaluation of HRCHY-CytoCommunity’s cross-platform generalization using the\nmouse hippocampus Slide-seq V2 and human breast cancer Visium V1 datasets.\n- Supplementary Figure 3. Spatial maps of fine-grained CNs in the human CRC Visium HD\ndataset.\n- Supplementary Figure 4. Spatial maps of coarse-grained TCs in the mouse hypothalamic\npreoptic region MERFISH dataset.\n- Supplementary Figure 5. Performance evaluation of HRCHY-CytoCommunity using the\nmouse ICH Stereo-seq dataset.\n- Supplementary Figure 6. Spatial maps of coarse-grained TCs in the mouse ICH Stereo-\nseq dataset.\n- Supplementary Figure 7. Cell-type map and cluster stability across different parameter\nsettings on the mouse hippocampus Slide-seq V2 dataset.\n- Supplementary Figure 8. Cell-type map and cluster stability across different parameter\nsettings on the human breast cancer Visium V1 dataset.\n- Supplementary Figure 11. Complementary results of survival analysis for the breast\ncancer IMC dataset.\n- Supplementary Figure 12. Average cell-type enrichment scores of unified CNs in the\nbreast cancer IMC dataset.\n- Supplementary Figure 13. Ablation study of key components in HRCHY-\nCytoCommunity using the simulated dataset.\n- Supplementary Figure 14. Ablation study of key components in HRCHY-\nCytoCommunity on coarse-grained TC identification using real spatial omics datasets.\n- Supplementary Figure 15. Ablation study of key components in HRCHY-\nCytoCommunity on fine-grained CN identification using real spatial omics datasets.\n- Supplementary Figure 16. Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity using the simulated dataset.\n- Supplementary Figure 17. Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity on coarse-grained TC identification using real spatial omics datasets.\n- Supplementary Figure 18. Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity on fine-grained CN identification using real spatial omics datasets.\n- Supplementary Figure 19. Robustness assessment of HRCHY-CytoCommunity to label\ninaccuracies and annotation resolution variations using the mouse hypothalamic preoptic\nregion MERFISH dataset.\n- Supplementary Figure 20. Robustness assessment of HRCHY-CytoCommunity to label\ninaccuracies and annotation resolution variations using the human TNBC MIBI-TOF\ndataset.\n- Supplementary Figure 21. Runtime benchmarking using Visium HD datasets.\n- Supplementary Figure 22. Memory usage benchmarking using Visium HD datasets.\n- Supplementary Figure 23. Performance evaluation of HRCHY-CytoCommunity using the\nsimulated dataset.\n- Supplementary Figure 24. Ambiguous structures in the mouse Slide-seq V2 dataset.\n- Supplementary Note 1. Fine-grained CNs identified in the mouse spleen CODEX and\nhuman CRC Visium HD datasets.\n- Supplementary Note 2. Coarse-grained TCs identified in the mouse hypothalamic preoptic\nregion MERFISH and intracerebral hemorrhage Stereo-seq datasets.\n- Supplementary Note 3. Ablation studies of HRCHY-CytoCommunity.\n- Supplementary Note 4. Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity.\n- Supplementary Note 5. Robustness assessment of HRCHY-CytoCommunity.\n- Supplementary Note 6. Scalability and speed analysis.\n- Supplementary Note 7. Experimental details for reproducibility.\n- Supplementary Note 8. Extension of HRCHY-CytoCommunity to three or more\nhierarchical levels.\n- Supplementary Note 9. Identification of ambiguous structures in hierarchical tissue\norganization.\n- Supplementary Table 1. Details of spatial omics datasets used for performance evaluation.\n- Supplementary Table 2. Hyperparameter settings of HRCHY-CytoCommunity used in all\ndatasets.\n- Supplementary Table 3. Hyperparameter settings of NeST used in all datasets.\n- Supplementary Table 4. Cell-type composition profiles for CNs in simulation data.\nIn the following, reviewers’ comments remain in black, our point-by-point answers appear in\nblue. Text and figures we have added or altered in the main manuscript are also highlighted in\nblue, while the supplementary information appears in black, as the entire document has been\ncomprehensively updated.\nResponse to reviewer comments\nReviewer #1 (Remarks to the Author):\nHRCHY-CytoCommunity presents a novel graph neural network (GNN)-based framework for\nidentifying hierarchical multicellular structures in spatial omics data. Leveraging differentiable\ngraph pooling and graph pruning, the method robustly uncovers nested tissue compartments\nand cellular neighborhoods from cell-type spatial maps. Overall, this work is a significant\nadvance in spatial omics analysis with biological and clinical insights.\nMajor advances of this approach:\n1) HRCHY-CytoCommunity outperforms the state-of-the-art NeST method across multiple\ndatasets (mouse spleen CODEX, hypothalamic MERFISH, TNBC MIBI-TOF). It accurately\nidentifies spatially discontinuous structures and achieves higher accuracy, adjusted mutual\ninformation, and Macro-F1 scores.\n2) The hierarchical majority voting strategy enhances reproducibility while preserving nested\nrelationships between TCs and CNs. Unlike NeST, it assigns all cells to structures, avoiding\ncoverage gaps.\n3) A cell-type enrichment-based clustering module aligns structures across samples, revealing\nconserved and patient-specific biological features (e.g., 26 unified CNs in TNBC).\n4) The method reveals functional hierarchies, such as neuron/non-neuron interactions in the\nhypothalamus and immune-tumor compartmentalization in TNBC, aligning with known\nbiology. Applied to breast cancer, HRCHY-CytoCommunity enables multi-scale prognostic\nstratification.\nLimitations:\nR1C1\n1. Researchers need to pre-define/select structure counts/levels. Are there any quantitative\nindexes to guide the hyper-parameter settings?\nWe appreciate the reviewer’s insightful comment. In the revised manuscript, inspired by the\nstate-of-the-art spatial domain detection method CellCharter1, we introduced a cluster stability\nmetric based on the Fowlkes-Mallows Index (FMI) to guide the selection of the numbers of\nfine-grained cellular neighborhoods (CNs, denoted as c₁) and coarse-grained tissue\ncompartments (TCs, denoted as c₂).\nIn this approach, users were first required to define a candidate search range for (c₁, c₂) and the\nnumber of repeated runs R. For each candidate pair (c₁, c₂), the model was trained R times\nindependently. The stability of the resulting hierarchical structures (i.e., CNs and TCs) was then\nquantified using the FMI computed across runs. Specifically, the FMI was computed separately\nfor the identified CNs and TCs across all pairs of runs for a given (c₁, c₂) configuration. This\nyielded a set of CN-level FMIs and a set of TC-level FMIs. For each level, the median FMI\nacross all run pairs was computed to summarize the cluster stability at that level. The final\nstability score for a (c , c ) configuration was defined as the average of the CN-level and TC-\n1 2\nlevel median FMIs. The optimal (c₁, c₂) pair was selected as the one that yields the highest\nstability score. Finally, among the R models trained under the optimal (c₁, c₂) configuration, the\nmodel achieving the lowest training loss was chosen to output the final hierarchical tissue\nstructures.\nThis procedure provides a quantitative criterion to guide the selection of (c₁, c₂) rather than\nrelying solely on ad hoc or domain-specific heuristics. More importantly, it also mitigates the\nrisk of over-segmentation when large values of (c₁, c₂) are used. Specifically, when the number\nof clusters is overestimated, the resulting clusters tend to be unstable across repeated runs,\nwhich is reflected by low FMI stability. Conversely, when the chosen (c₁, c₂) properly matches\nthe underlying tissue organization, the model yields reproducible partitions with high FMI\nstability.\nTaken together, instead of requiring users to specify (c₁, c₂) solely based on prior knowledge,\nHRCHY-CytoCommunity employs an FMI-based cluster stability criterion to determine the\noptimal values of (c₁, c₂), while guarding against over-segmentation. This ensures that HRCHY-\nCytoCommunity can identify robust and stable hierarchical structures.\nFull details are provided in the Methods section under the title “Determination of optimal\nnumbers of hierarchical tissue structures based on cluster stability”.\nR1C2\n2. Can this model be adapted to RNA-seq based spot-structure data, such as Visium HD &\nStereoSeq?\nWe appreciate the reviewer for raising this valuable point. Yes, HRCHY-CytoCommunity\ncan be directly applied to RNA-seq based spot-structure data such as Visium HD and\nStereo-seq. To demonstrate this, we systematically evaluated its performance on four\nrepresentative sequencing-based spatial transcriptomics datasets spanning diverse resolutions\nand technologies. Specifically:\n1. Human colorectal cancer (CRC) Visium HD dataset (8 µm-bin resolution): HRCHY-\nCytoCommunity was the only method that accurately distinguished tumor and normal\ncompartments in the P2CRC sample. Similarly, in the remaining samples, these two\ncompartments were also clearly delineated, though with moderate inclusion of intestinal\nepithelial-enriched regions within the tumor compartment (Response Fig. 1b, left column).\nQuantitatively, HRCHY-CytoCommunity achieved the highest average AMI and Macro-\nF1 scores among all representative state-of-the-art methods (Response Fig. 1d).\n2. Mouse intracerebral hemorrhage (ICH) Stereo-seq dataset (Single-cell resolution):\nHRCHY-CytoCommunity achieved high accuracy in identifying fine-grained CNs,\nparticularly excelling in detecting small-scale CNs such as olfactory areas (OLF), meninges,\nand pallidum (PAL) (Response Fig. 2b, left column). Quantitatively, HRCHY-\nCytoCommunity achieved significantly higher AMI and Macro-F1 scores than NeST (one-\nsided paired t-test, P-values < 1E-7) and other non-hierarchical methods (one-sided paired\nt-test, P-values < 0.001), with performance comparable to NicheCompass (Response Fig.\n2d).\n3. Mouse hippocampus Slide-seq V2 dataset (10 µm-spot resolution): HRCHY-\nCytoCommunity recovered biologically meaningful hierarchical structures. For coarse-\ngrained TCs, our model accurately delineated the cortical region, hippocampus, and brain\nstem (Response Fig. 3b, middle panel). For fine-grained CNs, our model accurately\nrecovered substructures including the CA1-CA3 regions, dentate gyrus, corpus callosum,\nand distinct cortical layers (Response Fig. 3c, middle panel), closely matching the reference\nAllen Mouse Brain Atlas2.\n4. Human breast cancer Visium V1 dataset (55 µm-spot resolution): HRCHY-\nCytoCommunity accurately separated tumor and non-tumor compartments (Response\nFig. 4b, middle panel) and identified fine-grained CNs consistent with manual annotation,\nsuch as DCIS/LCIS_1, DCIS/LCIS_3, and IDC_4 (Response Fig. 4c, middle panel), in line\nwith manual annotation from Xu et al.3.\nIn summary, these results demonstrate the strong generalizability of HRCHY-\nCytoCommunity to RNA-seq based spot-structure data spanning diverse technologies and\nresolutions.\nFull details of the extended performance evaluation on sequencing-based data are provided in\nthe Results section under the titles: “Performance evaluation on coarse-grained TC\nidentification”, “Performance evaluation on fine-grained CN identification”, and “Cross-\nplatform generalization of HRCHY-CytoCommunity to low-resolution spatial\ntranscriptomics”, as well as in Supplementary Notes 1 and 2.\nResponse Fig. 1 (extracted from Fig. 2). Performance evaluation of HRCHY-CytoCommunity on the\nhuman CRC Visium HD dataset. (a) Spatial maps of human CRC samples (8 μm bins) generated by Visium\nHD technology. Bins are colored by cell types (left) and manually annotated coarse-grained structures (right).\n(b) Coarse-grained TCs identified by HRCHY-CytoCommunity (left) and NeST (right). (c) Coarse-grained\nTCs identified by non-hierarchical methods, including CellCharter, NicheCompass, and SpaSEG. (d) AMI\nand Macro-F1 scores calculated using manual compartment annotations. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same sample are connected by grey dashed lines. P-values were calculated using\none-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. “Unidentified” indicates\ncells not assigned to any TC. “Unmatched” denotes TCs that could not be aligned with manual annotations.\nResponse Fig. 2 (Supplementary Fig. 5). Performance evaluation of HRCHY-CytoCommunity on the\nmouse ICH Stereo-seq dataset. (a) Spatial maps showing the distribution of cell types (left) and molecularly\ndefined regions (right). (b) Fine-grained CNs identified by HRCHY-CytoCommunity (left) and NeST (right).\n(c) Fine-grained CNs identified by non-hierarchical tissue structure identification methods, including\nCellCharter, NicheCompass, GraphST, and SpaGCN. “Unidentified” indicates cells not assigned to any CN.\n(d) AMI and Macro-F1 scores calculated using molecularly defined regions. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same samples are connected by grey dashed lines. P-values were calculated using\na one-sided paired t-test. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. “Unidentified” indicates\ncells not assigned to any CN. L2/3, cortical layer 2 and cortical layer 3; L4, cortical layer 4; L5, cortical layer 5;\nL6, cortical layer 6; Lesion, Lesion region; injured_STR, injured striatum region; injured_cortex, injured cortex\nregion; injured_PAL, injured pallidum region; DG, dentate gyrus; FT, fiber tracts; HY, hypothalamus; OLF,\nolfactory areas; MH, medial habenula; CA, Ammon’s horn; STRd, striatum dorsal region; STRv, striatum\nventral region; TH, thalamus; VS, ventricular systems; isl, islm, islands of Calleja and major island of Calleja;\nPAL, pallidum; RT, reticular nucleus of the thalamus; CAsp, field Ammon’s horn, pyramidal layer; LSX, lateral\nseptal complex; sAMY, striatum-like amygdala nuclei; CTXsp, cortical subplate.\nResponse Fig. 3 (extracted from Fig. 4 and Supplementary Fig. 6). Performance evaluation of HRCHY-\nCytoCommunity on the mouse hippocampus Slide-seq V2 dataset. (a) Cluster stability measure used by\nHRCHY-CytoCommunity to determine the optimal numbers of coarse-grained TCs and fine-grained CNs.\nThe Fowlkes-Mallows Index-based cluster stability values under each parameter setting are connected by\nsolid lines. The red circle marks the selected stable resolution of tissue structures. (b) Coarse-grained TC\nreference from the Allen Mouse Brain Atlas 2 (left), and TCs identified by HRCHY-CytoCommunity (middle)\nand NeST (right) in the mouse hippocampus. (c) Fine-grained CN reference from the Allen Mouse Brain\nAtlas (left), and CNs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the mouse\nhippocampus. V3, third ventricle; MH, medial habenula; LH, lateral habenula; LP, lateral posterior nucleus\nof the thalamus; LD, lateral dorsal nucleus of the thalamus. “Unidentified” indicates cells not assigned to any\nTC or CN. “Unmatched” denotes TCs or CNs that could not be aligned with the Allen Mouse Brain Atlas.\nResponse Fig. 4 (extracted from Fig. 4 and Supplementary Fig. 7). Performance evaluation of HRCHY-\nCytoCommunity on the human breast cancer Visium V1 dataset. (a) Cluster stability measure used by\nHRCHY-CytoCommunity to determine the optimal numbers of coarse-grained TCs and fine-grained CNs.\nThe Fowlkes-Mallows Index-based cluster stability values under each parameter setting are connected by\nsolid lines. The red circle marks the selected stable resolution of tissue structures. (b) Manual annotation of\ncoarse-grained TCs from Xu et al.3 (left), and TCs identified by HRCHY-CytoCommunity (middle) and\nNeST (right) in the human breast cancer. (c) Manual annotation of fine-grained CNs from Xu et al. (left),\nand CNs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the human breast cancer.\n“Unidentified” indicates cells not assigned to any TC or CN. Red arrowheads indicate the IDC_3 region\nmanually annotated by Xu et al., which was reclassified by both HRCHY-CytoCommunity and NeST as a\nnon-tumor CN.\nReviewer #2 (Remarks to the Author):\nIn this study, Wang et al. introduced a spatial clustering framework HRCHY-CytoCommunity.\nUnlike conventional spatial clustering methods such as SpaGCN or graphST that deliver a\nsingle level spatial domain clustering, the proposed method can present hierarchical spatial\nclustering, detecting nested spatial domains and structures that can not be detected with the\nconventional methods. Considering the complexity of the tissue microenvironment in various\ncomplex diseases such as cancer, the capability to detect hierarchical spatial clustering can help\nus better understand the spatial patterns underlying various biological conditions and diseases,\ntherefore it surely will be of interest of an ever-growing community of single-cell spatial\ngenomics. Although I am fully convinced the importance of the problem and appreciate the\nauthors effort in addressing this issue, I found that the authors may need to enhance the\nbenchmarking (current only one single method was compared and benchmarked) and add a\nablation study (to showcase the necessity of the all modules that they utilized). Beyond showing\nthe effectiveness of the method and its superior performance in comparison with other method,\nthe authors may also want to show more evidence to justify why the proposed touchup major\nvoting strategy (NOT a integral part of the model and thus can’t be optimized all together) is\nsolid and won’t damage the performance dramatically compared to other integrative alternative\nstrategy (see major comment #2)\nMajor concerns\nR2C1\n(1) A potential limitation of the HRCHY-CytoCommunity method lies in its reliance on user-\nspecified hyperparameters c1 and c2, which define the number of fine-grained cellular\nneighborhoods (CNs) and coarse-grained tissue compartments (TCs), respectively. While the\nauthors acknowledge that the model can learn fewer clusters than specified if overestimated,\nthere is no principled strategy provided for selecting appropriate values of c1 and c2 in practice.\nThis imposes a significant constraint on the method’s generalizability, particularly in real-world\nsettings where prior knowledge about hierarchical tissue architecture is limited or absent.\nWithout an automated or data-driven approach to determine optimal resolution, users must rely\non ad hoc or domain-specific heuristics, which may lead to suboptimal or biologically\nmisleading interpretations. A discussion on how to systematically tune or infer these\nhyperparameter, ideally integrated to be part of the model, would substantially improve the\napplicability and robustness of the framework.\nWe appreciate the reviewer’s insightful comment and acknowledge that requiring users to\npredefine the numbers of fine-grained cellular neighborhoods (CNs, denoted as c₁) and coarse-\ngrained tissue compartments (TCs, denoted as c₂) could constrain the generalizability of\nHRCHY-CytoCommunity. To address this issue, inspired by the state-of-the-art spatial domain\ndetection method CellCharter1, we introduced a data-driven heuristic strategy based on\ncluster stability to determine the optimal values of (c₁, c₂) in the revised manuscript.\nIn this approach, users were first required to define a candidate search range for (c₁, c₂) and the\nnumber of repeated runs R. For each candidate pair (c₁, c₂), the model was trained R times\nindependently. The stability of the resulting hierarchical structures (i.e., CNs and TCs) was then\nquantified using the Fowlkes-Mallows Index (FMI) computed across runs. Specifically, the\nFMI was computed separately for the identified CNs and TCs across all pairs of runs for a given\n(c₁, c₂) configuration. This yielded a set of CN-level FMIs and a set of TC-level FMIs. For each\nlevel, the median FMI across all run pairs was computed to summarize the cluster stability at that\nlevel. The final stability score for a (c , c ) configuration was defined as the average of the\n1 2\nCN-level and TC-level median FMIs. The optimal (c₁, c₂) pair was selected as the one that\nyields the highest stability score. Finally, among the R models trained under the optimal (c₁, c₂)\nconfiguration, the model achieving the lowest training loss was chosen to output the final\nhierarchical tissue structures.\nThis procedure provides a data-driven criterion to guide the selection of (c₁, c₂) rather than\nrelying solely on ad hoc or domain-specific heuristics. More importantly, it also mitigates the\nrisk of over-segmentation when large values of (c₁, c₂) are used. Specifically, when the number\nof clusters is overestimated, the resulting clusters tend to be unstable across repeated runs,\nwhich is reflected by lower FMI stability. Conversely, when the chosen (c₁, c₂) properly\nmatches the underlying tissue organization, the model yields reproducible partitions with high\nFMI stability.\nTaken together, instead of requiring users to specify (c₁, c₂) solely based on prior knowledge,\nHRCHY-CytoCommunity employs a data-driven heuristic strategy to determine the\noptimal values of (c₁, c₂), while guarding against over-segmentation. This ensures that\nHRCHY-CytoCommunity can identify robust and stable hierarchical tissue structures.\nFull details are provided in the Methods section under the title “Determination of optimal\nnumbers of hierarchical tissue structures based on cluster stability”\n.\nR2C2\n(2) The hierarchical majority voting strategy used in HRCHY-CytoCommunity, while helpful\nfor improving stability and enforcing a nested structure, is applied post hoc and is not integrated\ninto the model’s training objective. This disconnection introduces two key issues. First, the\nlearned soft assignments are not explicitly optimized to yield stable or hierarchically consistent\noutputs, which may lead to mismatches between training-time objectives and final outputs.\nSecond, because the voting procedure is non-differentiable, it prevents end-to-end training and\nmay obscure the interpretation of how hierarchical structures are learned. More critically, this\ntouch-up procedure could obscure systematic instability in the underlying GNN assignments,\nwhich would not be evident without the ensemble step. As a result, the biological or clinical\nsignificance of the identified structures may hinge on a potentially brittle post-processing routine.\nAs a practical alternative, the authors could consider incorporating a consistency regularization strategy, which encourages\nalignment of cluster assignments across stochastic runs during training, or a differentiable nesting regularization term, which\nsoftly penalizes misalignment between fine- and coarse-grained structures. These approaches would enable the model to\nlearn hierarchical relationships more robustly and transparently during optimization. One quick strategy is to introduce\nconsistency regularization during training. This involves running the model multiple times with stochastic perturbations (e.g.,\ndropout or graph augmentation) and penalizing discrepancies in the predicted cluster assignments across runs. By\nminimizing divergence (e.g., KL divergence) between outputs, the model is encouraged to produce stable and robust\nhierarchical assignments without relying on post hoc voting.\nAt a minimum, the authors should implement and compare one of these integrative alternatives\nagainst the current ensemble-based (touchup major voting based) strategy to demonstrate that\nthe model's performance and robustness are not overly dependent on post hoc correction. Such\na comparison would strengthen the claims of generalizability and biological interpretability. If\nthe performances are indeed similar (even slightly worse), then the authors may not need to\nupdate all the results with the new model structure, which is quite intensive.\nWe appreciate the reviewer for this constructive suggestion. We fully agree that the post hoc\nmajority-voting ensemble has limitations, as it is non-differentiable and not directly integrated\ninto the model training. In the revised manuscript, we have addressed this concern by replacing\nthe post hoc ensemble with a dropout-based consistency regularization strategy during\ntraining. Specifically,\n We introduced a consistency regularization term (adapted from the GRAND4 framework)\nthat enforces agreement of hierarchical structure assignments across multiple stochastic\nperturbations, where these perturbations are generated using DropNode (random masking\nof node attribute vectors ) (Response Fig. 5b).\n We further introduced an entropy-based balance regularization term that encourages\nbalanced cluster assignments by maximizing the entropy of the predicted distributions\nacross clusters to the prevent collapse of the coarse-grained TC assignments.\n These regularization terms were jointly optimized with the unsupervised hierarchical tissue\nstructure assignment objective, enabling fully end-to-end training without relying on\npost hoc corrections (Response Fig. 5b).\nResponse Fig. 5 (Fig. 1). Overview of the HRCHY-CytoCommunity framework. The framework\nformulates hierarchical tissue structure identification as a hierarchical community detection problem on a\ncell-cell proximity graph. HRCHY-CytoCommunity includes two modules: (a) a soft hierarchical structure\nassignment base module and (b) a robust hierarchical structure identification module. (a) The base module\nbegins with a cell-cell proximity graph, where nodes represent cells and their attributes are one-hot encoded\ncell-type vectors. A graph convolution layer followed by a fully-connected layer transforms node attribute\nvectors into soft fine-grained cellular neighborhood (CN) assignment vectors. A differentiable graph\npooling layer is then used to generate a coarsened completed graph, where each pooled node represents a\nCN. Adaptive edge pruning is applied to mitigate over-smoothing. Another graph convolution and fully-\nconnected layers process the pooled node embeddings to produce soft coarse-grained tissue compartment\n(TC) assignment vectors. The entire process is optimized by a graph minimum cut (MinCut)-based loss\nfunction, denoted as L . (b) To enhance the accuracy and robustness of hierarchical tissue structure\nBase\nidentification, M perturbed cell-cell proximity graphs are generated via DropNode (random masking of\nnode attribute vectors). Each graph is processed by the base module to produce soft fine-grained CN and\ncoarse-grained TC assignment matrices. Consistency and balance regularization during training enforces\nstability across these M assignments. The total loss, denoted as L , combines the base MinCut loss (L )\nTotal Base\nwith a consistency regularization term (L , which minimizes divergence across M assignments) and a\nConsis\nbalance regularization term (L , which prevents cluster collapse). Finally, this module yields stable\nBalance\nassignment matrices, from which final robust hierarchical tissue structures are derived.\nTo assess the effectiveness of the consistency and balance regularization terms, we conducted\nablation studies on both simulated and real datasets. The results demonstrated that HRCHY-\nCytoCommunity with these regularization terms substantially improves the accuracy and\nstability of coarse-grained TC identification compared to HRCHY-CytoCommunity\nwithout regularization (denoted as base model) (Response Figs. 6a and b). For fine-grained\nCN identification, HRCHY-CytoCommunity consistently achieved stable and strong\nperformance (Response Figs. 6a and c)\nTo further evaluate whether this regularization strategy provides additional advantages in\naccuracy over the original post hoc ensemble procedure, we compared the performance of the\nbase model trained with consistency and balance regularization (i.e., Base model + Consistency\n& Balance reg) against the majority voting ensemble of 10 independent runs of the base model\n(i.e., Base model + Majority voting ensemble). The results showed that the regularized model\nconsistently achieved higher accuracy across all real datasets (Response Fig. 7a).\nWe also examined whether the regularization strategy improves model stability and prevents\nthe cluster collapse by comparing the similarity of predictions across repeated independent runs\nin the base model and the regularized model. Stability was quantified using the FMI between\npairs of resulting clusters. The results showed that introducing the consistency and balance\nregularization terms enhances the model stability and effectively prevents cluster collapse\n(Response Figs. 7b and c).\nIn summary, these results demonstrate that incorporating the consistency and balance\nregularization terms substantially improves both accuracy and stability of HRCHY-\nCytoCommunity. Accordingly, in the revised manuscript, all ensemble-based results reported\nin the original version have been replaced by those obtained from the model with\nconsistency and balance regularization.\nFull details of the consistency and balance regularization strategy are provided in the\nMethods section under the title “Robust hierarchical tissue structure identification”. Full\ndetails of the ablation studies are provided in Supplementary Note 3 and Supplementary\nFigs. 13-15.\nResponse Fig. 6 (extracted from Supplementary Figs. 13-15). Ablation study of consistency and balance\nregularization on simulated and real spatial omics datasets. Hierarchical tissue structure identification\nunder different model configurations. Each configuration was run 10 times with different random seeds.\nBoxplot elements are defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5×\ninterquartile range. Performances are quantified by Adjusted Mutual Information (AMI) and Macro-F1\nscores. Each point corresponds to the performance of one run. Ablations were performed on (a) the simulated\ndataset, (b) real datasets with ground-truth coarse-grained TCs including the mouse spleen CODEX dataset\n(imaging-based) and the human colorectal cancer (CRC) Visium HD dataset (sequencing-based), as well as\n(c) real datasets with ground-truth fine-grained CNs including the mouse hypothalamic preoptic region\nMERFISH dataset (imaging-based) and the mouse intracerebral hemorrhage (ICH) Stereo-seq dataset\n(sequencing-based).\nResponse Fig. 7. Accuracy and stability evaluation of different versions of HRCHY-CytoCommunity.\n(a) Accuracy evaluation comparing the base model of HRCHY-CytoCommunity trained with consistency\nand balance regularization terms (i.e., Base model + Consistency & Balance reg) against the base model with\npost hoc ensemble (i.e., Base model + Majority voting ensemble) across all real datasets with ground-truth\nannotations generated from diverse technologies (i.e., the mouse spleen CODEX dataset, the human\ncolorectal cancer Visium HD dataset, the mouse hypothalamic preoptic region MERFISH dataset, and the\nmouse intracerebral hemorrhage Stereo-seq dataset). Accuracy is quantified by Adjusted Mutual Information\n(AMI) and Macro-F1 scores. (b-c) Stability evaluation of (b) identified coarse-grained TCs and (c) identified\nfine-grained CNs from the base model and the regularized model (i.e., Base model + Consistency & Balance\nreg).The stability score is quantified by Fowlkes-Mallows Index (FMI). Specifically, when all cells were\nassigned to a single cluster in a run (cluster collapse), the stability score involveing that run was set to zero.\nBoxplot elements are defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5×\ninterquartile range. Each dot represents the stability score between results from two independent runs with\nvarying random seeds.\nR2C3\n(3) The developed method integrates multiple components, including a hierarchical GNN with\ndifferentiable pooling, a post hoc majority-voting ensemble for robustness, and a clustering-\nbased module for cross-sample alignment. However, the manuscript does not provide an\nablation study to evaluate the individual contribution of these modules to overall performance.\nIn particular, it remains unclear how much of the reported improvement over baseline methods\narises from the GNN-based hierarchical learning itself, versus the ensemble voting or cross-\nsample unification steps. An ablation study removing or modifying each module, such as\ntraining without voting, disabling edge pruning, or skipping the alignment step, would be\nvaluable for confirming that each component is essential and for understanding where the gains\nin performance and robustness originate.\nWe appreciate the reviewer’s valuable suggestion. In the revised manuscript, we conducted\ncomprehensive ablation studies on both simulated and real datasets to evaluate the individual\ncontributions of key components in HRCHY-CytoCommunity. Notably, in this revised version,\nthe post hoc majority-voting ensemble used in the original manuscript has been replaced by the\nconsistency and balance regularization strategy, and the clustering-based cross-sample\nalignment module is an optional component that is not involved in performance benchmarking.\nTherefore, our ablation studies focus on the following three key components:\n1) Consistency and balance regularization terms, designed to improve the stability of\nhierarchical tissue structure assignments and prevent cluster collapse of coarse-grained TC\nidentification.\n2) Adaptive edge pruning, designed to automatically remove edges between fine-grained\nCNs to reduce over-smoothing during message passing in coarse-grained TC identification.\n3) Adaptive α scheduling, which dynamically balances fine- and coarse-level tissue structure\nidentification during training.\nOur results showed that:\n1. For the task of coarse-grained TC identification:\na) The best performance was consistently achieved when both the consistency and\nbalance regularization terms were applied together, confirming their complementary\nroles in improving accuracy and stability while preventing cluster collapse (Response\nFigs. 8a and 9a).\nb) The adaptive edge-pruning strategy outperformed both extreme configurations (i.e.,\npreserving all edges and removing all edges) in both simulated and real datasets,\nconfirming its effectiveness in filtering trivial connections between fine-grained CNs\n(Response Figs. 8b and 9b).\nc) The adaptive α scheduling strategy had little effect on the simulated dataset but\nimproved performance in real datasets, suggesting that it facilitates progressive\noptimization of hierarchical clustering under complex biological structures (Response\nFigs. 8c and 9c).\n2. For the task of fine-grained CN identification:\na) The three components did not substantially alter performance, which remained stable\nacross different settings in both simulated and real datasets (Response Figs. 8 and 10).\nThis is consistent with expectations, as these components were primarily designed for\npreventing suboptimal coarse-grained TC identification.\nTaken together, these ablation studies demonstrate that the three evaluated components of\nHRCHY-CytoCommunity play critical and complementary roles in coarse-grained TC\nidentification. In contrast, for fine-grained CN identification, the model maintains robust\nperformance across varying component settings.\nFull details of the ablation studies are provided in Supplementary Note 3 and Supplementary\nFigs. 13-15.\nResponse Fig. 8 (Supplementary Fig. 13). Ablation study of key components in HRCHY-\nCytoCommunity using the simulated dataset. Hierarchical tissue structure identification under different\nmodel configurations. Each configuration was run 10 times with different random seeds. Boxplot elements\nare defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range.\nPerformances are quantified by Adjusted Mutual Information (AMI) and Macro-F1 scores. Each point\ncorresponds to the performance of one run. Ablation studies were performed on (a) consistency and balance\nregularization terms, (b) the adaptive edge-pruning strategy, and (c) the adaptive α scheduling strategy.\nResponse Fig. 9 (Supplementary Fig. 14). Ablation study of key components in HRCHY-\nCytoCommunity on coarse-grained TC identification using real spatial omics datasets. Hierarchical\ntissue structure identification under different model configurations using the mouse spleen CODEX dataset\n(imaging-based) and the human colorectal cancer (CRC) Visium HD dataset (sequencing-based). Each run\nwas run 10 times with different random seeds. Boxplot elements are defined as: center line, median; box\nlimits, upper and lower quartiles; whiskers, 1.5× interquartile range. Performances are quantified by\nAdjusted Mutual Information (AMI) and Macro-F1 scores. Each point corresponds to the performance of\none run. Ablation studies were performed on (a) consistency and balance regularization terms, (b) the\nadaptive edge-pruning strategy, and (c) the adaptive α scheduling strategy.\nResponse Fig. 10 (Supplementary Fig. 15). Ablation study of key components in HRCHY-\nCytoCommunity on fine-grained CN identification using real spatial omics datasets. Hierarchical tissue\nstructure identification under different model configurations using the mouse hypothalamic preoptic region\nMERFISH dataset (imaging-based) and the mouse intracerebral hemorrhage (ICH) Stereo-seq dataset\n(sequencing-based). Each configuration was run 10 times with different random seeds. Boxplot elements\nare defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range.\nPerformances are quantified by Adjusted Mutual Information (AMI) and Macro-F1 scores. Each point\ncorresponds to the performance of one run. Ablation studies were performed on (a) consistency and balance\nregularization terms, (b) the adaptive edge-pruning strategy, and (c) the adaptive α scheduling strategy.\nR2C4\n(4) The current benchmarking is restricted to a single method, NeST, which, while hierarchical,\nis designed for high-dimensional spatial transcriptomics and does not generalize well to low-\ndimensional or cell-type-based spatial data such as spatial proteomics. This narrow comparison\nlimits the ability to assess the broader utility and competitiveness of HRCHY-CytoCommunity.\nIn particular, several more recent and widely used frameworks—such as SpaSEG, Giotto, and\nGraphST—are highly relevant comparators. Although these methods do not explicitly model\nhierarchical structures, they can be run at multiple resolutions or applied recursively to simulate\nhierarchical segmentation. Benchmarking against such methods would provide a more\ncomprehensive evaluation of whether HRCHY-CytoCommunity's explicit hierarchical\nmodeling offers practical advantages in terms of spatial coherence, robustness, and biological\nrelevance. The authors are strongly encouraged to include comparisons to these additional\nmethods or, at a minimum, to discuss their relative performance.\nWe thank the reviewer for the important suggestion. Following the suggestion from Reviewers\n2 and 3, in the revised manuscript, we have benchmarked HRCHY-CytoCommunity against six\nstate-of-the-art non-hierarchical tissue structure identification methods, including SpaSEG 5,\nGiotto Suite 6, GraphST 7, SpaGCN 8, NicheCompass 9, and CellCharter 1. Since these\nmethods do not explicitly identify hierarchical tissue structures, we applied them separately at\ncoarse-grained TC and fine-grained CN resolutions. Performances were evaluated using\nAdjusted Mutual Information (AMI) and Macro-F1 scores across two datasets with ground-\ntruth TCs and two datasets with ground-truth CNs.\nOur results demonstrate that HRCHY-CytoCommunity consistently achieved superior or\ncompetitive performance than non-hierarchical methods. Specifically:\n1. For coarse-grained TC identification:\na) On the mouse spleen CODEX dataset (imaging-based), HRCHY-CytoCommunity\nwas the only method that consistently resolved the biologically meaningful lymphoid\ncompartment containing B-zone, periarteriolar lymphoid sheath, and marginal zone\nacross all samples (Response Fig. 11b, left column), whereas non-hierarchical methods\nlargely failed to capture this hierarchical structure (GraphST succeeded only in one\nsample). Quantitatively, HRCHY-CytoCommunity achieved significantly higher AMI\nand Macro-F1 scores than most non-hierarchical methods (one-sided paired t-test P-\nvalues < 0.05), while exhibiting performance comparable to GraphST (Response Fig.\n11d).\nb) On the human colorectal cancer (CRC) Visium HD dataset (sequencing-based),\nHRCHY-CytoCommunity was the only method that accurately distinguished tumor\nand normal compartments in the P2CRC sample, while all non-hierarchical methods\nfailed to do so. Similarly, in the remaining samples, these two compartments were also\nclearly delineated, though with moderate inclusion of intestinal epithelial-enriched\nregions within the tumor compartment (Response Fig. 12b, left column).\nQuantitatively, HRCHY-CytoCommunity achieved the highest average AMI and\nMacro-F1 scores among all representative non-hierarchical methods (Response Fig.\n12d).\n2. For fine-grained CN identification:\na) On the mouse hypothalamic preoptic region MERFISH dataset (imaging-based),\nHRCHY-CytoCommunity successfully recovered CNs consistent with manually\ndelineated nuclei with clear boundaries (Response Fig. 13b, left column), whereas\nother non-hierarchical methods identified intermixed CNs and consistently failed to\ncapture smaller regions like PS and PaAP in most samples (Response Fig. 13c).\nQuantitatively, HRCHY-CytoCommunity significantly outperformed all non-\nhierarchical methods (one-sided paired t-test, p < 0.05) in both AMI and Macro-F1\nmetrics, with the exception of CellCharter on the AMI and GraphST on the Macro-F1\nscore, where the differences were not statistically significant (Response Fig. 13d).\nb) On the mouse intracerebral hemorrhage (ICH) Stereo-seq dataset (sequencing-\nbased) HRCHY-CytoCommunity achieved high accuracy in identifying fine-grained\nCNs, particularly excelling in detecting small-scale CNs such as olfactory areas (OLF),\nmeninges, and pallidum (PAL) (Response Fig. 14b, left column). Quantitatively,\nHRCHY-CytoCommunity achieved significantly higher AMI and Macro-F1 scores\nthan most non-hierarchical methods (one-sided paired t-test, P-values < 0.001), with\nperformance comparable to NicheCompass (Response Fig. 14d).\nIn summary, these results indicate that HRCHY-CytoCommunity is not only quantitatively\ncompetitive but also qualitatively superior, as it identifies biologically meaningful structures\nmissed by non-hierarchical models. This highlights that the explicit hierarchical framework\nyields more accurate, stable, and interpretable multi-scale tissue organization across\ndiverse technologies.\nFull implementation details of these non-hierarchical spatial domain detection methods are\nprovided in the Methods section under the title “Running of published methods”. Full details\nof the performance evaluation are provided in the Results section under the titles:\n“Performance evaluation on coarse-grained TC identification”, and “Performance\nevaluation on fine-grained CN identification”, as well as in Supplementary Notes 1 and 2.\nResponse Fig. 11 (extracted from Fig. 2) Performance evaluation of HRCHY-CytoCommunity on the\nmouse spleen CODEX dataset. (a) Spatial maps of healthy mouse spleen samples generated by CODEX\ntechnology. Cells are colored by cell types (left) and manually annotated compartments (right), which serve\nas the ground truth for coarse-grained TCs. (b) Coarse-grained TCs identified by HRCHY-CytoCommunity\n(left) and NeST (right). (c) Coarse-grained TCs identified by non-hierarchical methods, including CellCharter,\nNicheCompass, GraphST, and SpaGCN. (d) Adjusted Mutual Information (AMI) and Macro-F1 scores\ncalculated using manual compartment annotations. Each point corresponds to the performance on an\nindividual sample, with horizontal bars indicating the mean performance across all samples. Points from the\nsame sample are connected by grey dashed lines. P-values were calculated using one-sided paired t-tests. *,\nP-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. “Unidentified” indicates cells not assigned to any\nTC. “Unmatched” denotes TCs that could not be aligned with manual annotations.\nResponse Fig. 12 (extracted from Fig. 2). Performance evaluation of HRCHY-CytoCommunity on the\nhuman CRC Visium HD dataset. (a) Spatial maps of human CRC samples (8 μm bins) generated by Visium\nHD technology. Bins are colored by cell types (left) and manually annotated coarse-grained structures (right).\n(b) Coarse-grained TCs identified by HRCHY-CytoCommunity (left) and NeST (right). (c) Coarse-grained\nTCs identified by non-hierarchical methods, including CellCharter, NicheCompass, and SpaSEG. (d) AMI\nand Macro-F1 scores calculated using manual compartment annotations. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same sample are connected by grey dashed lines. P-values were calculated using\none-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. “Unidentified” indicates\ncells not assigned to any TC. “Unmatched” denotes TCs that could not be aligned with manual annotations.\nResponse Fig. 13 (Fig. 3). Performance evaluation of HRCHY-CytoCommunity on the mouse\nhypothalamic preoptic region MERFISH dataset. (a) Spatial maps of mouse hypothalamic preoptic\nregion samples generated by MERFISH technology. Cells are colored by cell types (left) and manually\nannotated hypothalamic nuclei (right), which were generated based on the outlines of hypothalamic nuclei\n(middle) and serve as the ground truth for fine-grained CNs. (b) Fine-grained CNs identified by HRCHY-\nCytoCommunity (left) and NeST (right). (c) Fine-grained CNs identified by non-hierarchical methods,\nincluding CellCharter, NicheCompass, GraphST, SpaGCN, SpaSEG, and Giotto Suite. (d) Adjusted Mutual\nInformation (AMI) and Macro-F1 scores calculated using manually annotated hypothalamic nuclei. Each\npoint corresponds to the performance on an individual sample, with horizontal bars indicating the mean\nperformance across all samples. Points from the same sample are connected by grey dashed lines. P-values\nwere calculated using one-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.\n“Unidentified” indicates cells not assigned to any CN.\nResponse Fig. 14 (Supplementary Fig. 5). Performance evaluation of HRCHY-CytoCommunity on the\nmouse ICH Stereo-seq dataset. (a) Spatial maps showing the distribution of cell types (left) and molecularly\ndefined regions (right). (b) Fine-grained CNs identified by HRCHY-CytoCommunity (left) and NeST (right).\n(c) Fine-grained CNs identified by non-hierarchical tissue structure identification methods, including\nCellCharter, NicheCompass, GraphST, and SpaGCN. “Unidentified” indicates cells not assigned to any CN.\n(d) AMI and Macro-F1 scores calculated using molecularly defined regions. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same samples are connected by grey dashed lines. P-values were calculated using\na one-sided paired t-test. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. L2/3, cortical layer 2\nand cortical layer 3; L4, cortical layer 4; L5, cortical layer 5; L6, cortical layer 6; Lesion, Lesion region;\ninjured_STR, injured striatum region; injured_cortex, injured cortex region; injured_PAL, injured pallidum\nregion; DG, dentate gyrus; FT, fiber tracts; HY, hypothalamus; OLF, olfactory areas; MH, medial habenula;\nCA, Ammon’s horn; STRd, striatum dorsal region; STRv, striatum ventral region; TH, thalamus; VS,\nventricular systems; isl, islm, islands of Calleja and major island of Calleja; PAL, pallidum; RT, reticular\nnucleus of the thalamus; CAsp, field Ammon’s horn, pyramidal layer; LSX, lateral septal complex; sAMY,\nstriatum-like amygdala nuclei; CTXsp, cortical subplate.\nR2C5\n(5) The reliance on cell-type annotations as input features introduces a dependence on upstream\npreprocessing quality, which is not evaluated. The authors should assess how robust their\nmethod is to label inaccuracies or varying annotation resolutions, or clarify the assumptions\nabout cell-type consistency across samples.\nWe thank the reviewer for raising this critical point. In the revised manuscript, we\nsystematically evaluated the robustness of HRCHY-CytoCommunity to both label inaccuracies\nand varying annotation resolutions using two datasets with complex cell-type compositions (i.e.,\nthe mouse hypothalamic preoptic region MERFISH dataset and the human TNBC MIBI-TOF\ndataset). Specifically:\n1. Robustness to label inaccuracies: We used the hierarchical tissue structures identified\nfrom unperturbed cell-type annotations as the reference. Cell-type labels were then\nrandomly perturbed by 10%–50%, and the accuracy and consistency of the resulting\nhierarchical structures were compared with the reference. The FMI metric was used to\nquantify the consistency between the perturbed and reference structures across different\nnoise levels.\na) In the MERFISH dataset, the identified hierarchical structures remained highly\nconsistent until the label error rate exceeded 40% (Response Fig. 15a; average FMI >\n0.5 when label error rate = 40%). The accuracy of the fine-grained CNs (measured\nby AMI and Macro-F1 scores) showed only a moderate decline, with performance\nlargely preserved under reasonable noise levels (Response Figs. 15b and c)\nb) In the MIBI-TOF dataset, hierarchical structures remained highly consistent across\nperturbation rates (Response Fig. 16a; average FMI > 0.65 when label error rate =\n50%). The accuracy of the identified coarse-grained TCs (measured by the proportion\nof neoplastic and immune cells correctly assigned to their corresponding TCs) was\nlargely preserved under reasonable levels of annotation noise, with only a moderate\ndecrease as the label error rate increased (Response Fig. 16b).\n2. Robustness to varying annotation resolutions: We further assessed robustness to\nannotation resolution by comparing resulting structures obtained under different levels of\ncell-type granularity (i.e., 15 vs. 9 cell types in MERFISH, 17 vs. 11 cell types in TNBC\nMIBI-TOF). The FMI was used to quantify the similarity between identified structures.\nSpecifically,\na) HRCHY-CytoCommunity consistently recovered stable hierarchical tissue structures\nacross annotation resolutions, achieving average FMI scores of 0.673 in the\nMERFISH dataset (Response Figs. 15d-f) and 0.590 in the MIBI-TOF dataset across\nsamples (Response Figs. 16c-e).\nTaken together, these analyses demonstrate that HRCHY-CytoCommunity is strongly robust\nto both label inaccuracies and variations in annotation resolution, alleviating concerns\nregarding the strong dependence on upstream cell type annotation quality.\nFull details of the robustness assessment are provided in Supplementary Note 5 and\nSupplementary Figs. 19 and 20.\nResponse Fig. 15 (Supplementary Fig. 19). Robustness assessment of HRCHY-CytoCommunity to\nlabel inaccuracies and annotation resolution variations using the mouse hypothalamic preoptic region\nMERFISH dataset. (a) Hierarchical tissue structure similarity under increasing cell-type label perturbation\nrates (10%-50%), quantified by Fowlkes-Mallows Index (FMI)-based robustness score. Higher values\nindicate greater structural consistency with the unperturbed reference. (b-c) Fine-grained CN identification\naccuracy under label perturbations, measured by Adjusted Mutual Information (AMI) (b) and Macro-F1\nscores (c). Each point represents performance on one sample, with horizontal bars represent the mean\nperformance across n = 5 samples. Performances (points) on the same sample are connected by grey dashed\nlines. Mean values under each label error rate are connected by blue solid lines. Shaded areas denote standard\ndeviation. (d) Single-cell spatial maps showing the distribution of 15 cell types (left), coarse-grained TCs\n(middle) and fine-grained CNs (right) identified by HRCHY-CytoCommunity based on this high-resolution\nannotation. (e) Single-cell spatial maps showing the distribution of 9 cell types (left), coarse-grained TCs\n(middle) and fine-grained CNs (right) identified by HRCHY-CytoCommunity based on this low-resolution\nannotation. (f) Hierarchical tissue structure similarity between different cell-type annotation resolutions,\nquantified by FMI-based robustness score. The red dashed line indicates the mean similarity (0.673) relative\nto the 15-cell-type annotation-based reference.\nResponse Fig. 16 (Supplementary Fig. 20). Robustness assessment of HRCHY-CytoCommunity to\nlabel inaccuracies and annotation resolution variations using the human TNBC MIBI-TOF dataset.\n(a) Hierarchical tissue structure similarity under increasing cell-type label perturbation rates (10%-50%),\nquantified by Fowlkes-Mallows Index (FMI)-based robustness score. Higher values indicate greater\nstructural consistency with the unperturbed reference. (b) Coarse-grained TC identification accuracy under\nlabel perturbations, measured by proportion of neoplastic and immune cells correctly assigned to their\ncorresponding TCs. Each point represents performance on one compartmentalized tumor sample, with\nhorizontal bars represent the mean performance across n = 15 samples. Performances (points) on the same\nsample are connected by grey dashed lines. Mean values under each label error rate are connected by blue\nsolid lines. Shaded areas denote standard deviation. (c) Single-cell spatial maps showing the distribution of\n17 cell types (left), coarse-grained TCs (middle) and fine-grained CNs (right) identified by HRCHY-\nCytoCommunity based on this high-resolution annotation. (d) Single-cell spatial maps showing the\ndistribution of 11 cell types (left), coarse-grained TCs (middle) and fine-grained CNs (right) identified by\nHRCHY-CytoCommunity based on this low-resolution annotation. (e) Hierarchical tissue structure\nsimilarity between different cell-type annotation resolutions, quantified by FMI-based robustness score.\nThe red dashed line indicates the mean similarity (0.590) relative to the 17-cell-type annotation-based\nreference.\nR2C6\n(6) Given the multi-step pipeline and ensemble-based framework, the scalability of HRCHY-\nCytoCommunity in terms of runtime and memory is not discussed. The authors should provide\nruntime/memory profiling or clarify the computational requirements to facilitate practical\nadoption.\nWe thank the reviewer for this valuable suggestion. In the revised manuscript, we have added\na dedicated scalability analysis to assess the runtime and memory efficiency of HRCHY-\nCytoCommunity. Specifically, we downsampled the largest Visium HD sample (P2CRC;\n473,318 spots, 18,085 genes) to generate subsets of 1k, 10k, 100k, and 473k spots. HRCHY-\nCytoCommunity was evaluated in two configurations: (i) the base model (without consistency\nand balance regularization), and (ii) the full model (with both regularization terms). We also\ncompared its performance against six other tissue structure identification methods, including\nNeST, CellCharter, NicheCompass, GraphST, SpaGCN, and SpaSEG. Our results showed that:\n1. Runtime scalability: Both models of HRCHY-CytoCommunity exhibited excellent\ncomputational efficiency. The base model achieved the fastest training speed, completing\ntraining and hierarchical structure identification on the full 473k-spot dataset in 1.39\nminutes. The full model, which involves repeated stochastic perturbations for stability,\ncompleted the same dataset in 8.02 minutes, ranking third after the base model and SpaSEG\n(Response Fig. 17a).\n2. Memory scalability: Both models of HRCHY-CytoCommunity showed nearly linear\ngrowth in memory usage with data size increasing, reflecting great scalability. On the full\n473k-spot dataset, the base model required only 4.09 GB of total memory, and the full\nmodel required 8.69 GB (Response Fig. 17b). A detailed breakdown further showed that\nGPU memory usage remained moderate for both models (2.70 GB for the base model and\n7.26 GB for the full model at 473k spots), indicating that HRCHY-CytoCommunity can be\nefficiently trained on commonly available GPUs, even for large-scale datasets (Response\nFig. 17d).\nIn summary, these results demonstrate that HRCHY-CytoCommunity achieves highly\ncompetitive scalability in both runtime and memory efficiency, supporting its practical\napplication to increasingly large and high-resolution spatial omics datasets.\nFull details of the scalability analysis are provided in Supplementary Note 6 and\nSupplementary Figs. 21 and 22.\nResponse Fig. 17 (extracted from Supplementary Figs 21 and 22) Runtime and memory scalability\nanalysis using the Visium HD dataset (a-b) Total runtime (a) and total memory usage (b) of eight methods\nincluding HRCHY-CytoCommunity without consistency and balance regularization (denoted as base model),\nHRCHY-CytoCommunity with consistency and balance regularization (denoted as full model), NeST,\nCellCharter, NicheCompass, GraphST, SpaGCN, and SpaSEG across Visium HD datasets of increasing sizes\n(from 1,000 to 473,318 spots) (c-d) Detailed decomposition of (c) runtime (d) memory usage for the all\nevaluated methods on the 473k-spot dataset. For each method, the runtime is separated into the training stage\n(red) and clustering stage (yellow); the total memory usage is separated into the GPU memory usage (green)\nand RAM usage (brown).\nMinor comments\nR2C1\n(1) HRCHY-CytoCommunity models a fixed two-level hierarchy. This assumption may limit\napplicability to tissues with deeper or non-uniform hierarchical organization. The authors\nshould discuss whether and how the framework can be extended to more flexible or adaptive\nhierarchical depths.\nWe thank the reviewer for this insightful comment. Indeed, the current implementation of\nHRCHY-CytoCommunity models a fixed two-level hierarchy (fine-grained CNs and coarse-\ngrained TCs), which may limit its applicability to tissues with deeper or non-uniform\nhierarchical organization. Conceptually, the framework is not restricted to two levels. By\niteratively applying the differentiable graph pooling module and extending the MinCut-based\nloss function, HRCHY-CytoCommunity can be naturally generalized to three or more\nhierarchical levels (e.g., smaller-scale niche–CN–TC). In this study, we focused on a two-level\nhierarchy due to the lack of reliable multi-level ground-truth references for quantitative\nbenchmarking, but the framework itself readily supports deeper hierarchies. For the non-\nuniform hierarchical organization, we agree that adaptive identification of tissue structures with\nvariable hierarchical depths represents an important direction for future development. In the\nrevised manuscript, we provided a detailed unified framework of extending HRCHY-\nCytoCommunity to arbitrary hierarchical depths in Supplementary Note 8 (Lines 285-331,\nPages 8-9 in Supplementary Information). We have also explicitly acknowledged this\nlimitation and outlined future directions in the revised Discussion section.\n“Lastly, although this study focused on a two-level hierarchy, the underlying framework of\nHRCHY-CytoCommunity can be readily extended to deeper hierarchies (Supplementary Note\n8). Future work may explore adaptive mechanisms to automatically infer the optimal\nhierarchical depth or to model tissues with non-uniform organizational principles”\nR2C2\n(2) The interpretability of the discovered tissue structures is not systematically evaluated\nbeyond cell-type enrichment. It would be valuable to assess how consistently the model\ndiscovers known anatomical or functional regions across samples and tissues, and how\nuninterpretable clusters are handled.\nWe thank the reviewer for this constructive suggestion. In the revised manuscript, we have\nsystematically evaluated the interpretability and consistency of hierarchical tissue structures\nidentified by HRCHY-CytoCommunity across multiple datasets with known anatomical or\nfunctional organization.\nFor tissues with well-characterized hierarchical architectures, including the mouse spleen\n(CODEX), mouse hypothalamus (MERFISH), mouse hippocampus (Slide-seq V2), and mouse\nbrain (Stereo-seq) datasets, HRCHY-CytoCommunity consistently identified anatomically and\nfunctionally meaningful regions across different samples. For example, in the spleen, two major\ncompartments corresponding to the red pulp and lymphoid compartments were robustly\nrecovered in all three samples. In the hypothalamus, the model consistently detected symmetric\nnuclei corresponding to distinct functional domains across five samples. Similarly, in the\nhippocampus dataset, HRCHY-CytoCommunity accurately delineated the cortical region,\nhippocampus, brain stem, and their corresponding substructures, in agreement with the known\nanatomical organization described in the Allen Mouse Brain Atlas 2. Moreover, for\ncompartmentalized tumor tissues, including colorectal cancer (Visium HD) and breast cancer\n(Visium V1 and MIBI-TOF), HRCHY-CytoCommunity consistently identified the tumor-\ndominated and normal/immune cell-dominated compartments.\nFor tissues lacking well-established histological ground truth, we annotated each identified TC\nor CN using cell-type enrichment analysis to determine the dominant cell populations and their\nassociated biological functions. This approach enables functional interpretation of previously\nuncharacterized clusters.\nCollectively, these analyses demonstrate that HRCHY-CytoCommunity consistently\nidentifies biologically interpretable and reproducible hierarchical tissue structures across\ndiverse tissue types and experimental platforms.\nR2C3\n(3) a 2-levels hierarchical structure”, “a 2-level hierarchical structure”\nWe thank the reviewer for this correc⇒tion. In the revised manuscript, we have revised the text\nfrom \"a 2-levels hierarchical structure\" to \"a 2-level hierarchical structure\" throughout the\nmanuscript.\nR2C4\n(4) Use of inconsistent capitalization for “Cellular Neighborhoods” and “Tissue Compartments”\n(sometimes capitalized, sometimes not). Maybe the authors can pick one consistent style (e.g.,\nlowercase unless at the start of a sentence).\nWe thank the reviewer for pointing out this inconsistency. In the revised manuscript, we have\nstandardized the terminology such that the full terms \"cellular neighborhoods\" and \"tissue\ncompartments\" are written in lowercase throughout the text, except when appearing at the\nbeginning of a sentence or in section titles. Their abbreviations (CNs and TCs) remain\ncapitalized, consistent with common conventions. This revision has been applied consistently\nacross the entire manuscript.\nR2C5\n(5) Line 179, “the coarse-grained red pulp compartment could be further divided into two fine-\ngrained substructures, red pulp and marginal zone”. This phrasing is somewhat confusing, as it\nsuggests that “red pulp” is both a coarse-grained and a fine-grained region, and it categorizes\nthe marginal zone, which is typically considered a transitional or fine-grained region between\nred and white pulp, as a substructure of red pulp. This may need clarification or correction,\nespecially for readers with histological expertise.\nWe thank the reviewer for bringing this up and apologize for the confusion. In the revised\nmanuscript, we have clarified the description and corrected the compartment terminology to\nred pulp compartment and lymphoid compartment (comprising the periarteriolar lymphoid\nsheath, B-zone, and marginal zone), in alignment with previous histological studies10–12. Some\nstudies collectively classify these three substructures as part of the white pulp, whereas others\ndescribe the marginal zone as a transitional structure between red and white pulp. To avoid\nambiguity and ensure biological accuracy, we therefore refer to this compartment as the\nlymphoid compartment throughout the revised manuscript.\nReviewer #2 (Remarks on code availability):\nThe github repository generally looks good to me.\nThe provided GitHub repository includes a basic example dataset (mouse brain MERFISH) and\na structured walkthrough consisting of four steps: graph construction, soft clustering, ensemble\nvoting, and visualization. While this demonstrates the core workflow, there are several issues.\nThe example dataset is limited in scope, benchmark datasets from the paper (e.g., spleen or\nTNBC) are not provided. Also, the documentation lacks clarity on hyperparameters, and the\nscripts have limited error handling.\nWe thank the reviewer for this valuable feedback on our code. In the revised version, we have\nsubstantially improved the GitHub repository to enhance the usability and reproducibility of\nHRCHY-CytoCommunity.\nSpecifically, we have made the following updates:\n1. Expanded the example data and tutorials. In addition to the basic mouse brain\nMERFISH dataset, we now provide two additional representative subsets of the mouse\nspleen (CODEX) and the human TNBC (MIBI-TOF) datasets used in the paper, along\nwith corresponding preprocessed inputs for testing the hierarchical structure\nidentification workflow. The full datasets can be accessed through the data portals\nreferenced in the Methods section, following their respective data access policies.\n2. Improved documentation. We have added detailed explanations of key\nhyperparameters, including DropNode rate, λ_consis, λ_balance, and number of\nperturbations, together with recommended values and tuning strategies.\n3. Enhanced robustness of scripts. Error handling and input validation have been\nincorporated to ensure that the code runs smoothly across different computational\nenvironments.\nThese updates ensure that users can reproduce the full HRCHY-CytoCommunity workflow—\nfrom graph construction to hierarchical structure identification—using publicly available data\nand well-documented parameter settings. To further enhance usability, we have packaged\nHRCHY-CytoCommunity as a Python library and released it on PyPI\n(https://pypi.org/project/hrchy-cytocommunity/). The updated source code for HRCHY-\nCytoCommunity is available at https://github.com/huBioinfo/HRCHY-CytoCommunity.\nReviewer #3 (Remarks to the Author)：\nThe manuscript reviewed introduces HRCHY-CytoCommunity, a method for identifying\nhierarchical tissue structures using only spatial graphs and cell-type annotations. The study\npresents a method to perform hierarchical structure identification using graph neural networks\nwith pruning and pooling. The study exemplified this strategy in multiple datasets, highlighting\nthe identification of discontinuous spatial patterns and presenting an improved performance in\ndatasets with few expression features (e.g., protein-based data).\nWhile these are promising goals, the manuscript relies heavily on the assumption that\nhierarchical organization in tissues can be fully captured using cell-type annotations alone. This\nmay be an oversimplification of complex tissues where cellular behavior is driven by signaling\ngradients, cellular interactions, or dynamic gene programs. We believe this together with others\nexposed below, limit the current version of the manuscript from being suitable for publication.\nIn our opinion, only a deeply revised version of the manuscript, with clarification of the\nfollowing concerns should be considered for publication.\nMajor Comments\nR3C1\n1. Assumption of Hierarchy: The method assumes that tissue structures are organized\nhierarchically, which is not universally true. While this may apply to some biological systems\n(e.g., lymphoid tissues, brain), in others, spatial patterns are more continuous or transitional.\nThe manuscript should explicitly discuss the biological validity of this assumption and its\nlimitations.\nWe appreciate the reviewer for this constructive suggestion. We fully agree that the assumption\nof hierarchical organization does not universally apply to all tissues. Indeed, our framework is\nmotivated by the observation that many biological systems are naturally organized in a\nhierarchical manner 10,13–15. For example, the spleen comprises coarse-grained compartments\n(i.e., red pulp and lymphoid pulp) that further contain fine-grained substructures (i.e., marginal\nzone, B-cell zone, and periarteriolar lymphoid sheath). Similarly, tumor ecosystems often\nexhibit multi-scale organization, from macro-tissue compartments 15–18 (e.g., tumor core and\nstroma) to micro-cellular neighborhoods 19–22 (e.g., tertiary lymphoid structures and\nmicropapillary patterns). Capturing such hierarchies facilitated the association between spatial\ncontext and tumor development, prognosis, and therapy response 15,23,24.\nFollowing the reviewer’s suggestion, we have clarified the scope of our method in the revised\nmanuscript. Specifically, in the Introduction section, we highlight several organs known to\nexhibit hierarchical organization, such as the brain, tumor ecosystems, and lymphoid tissues\nlike the spleen, to better contextualize the motivation for HRCHY-CytoCommunity. In the\nDiscussion section, we explicitly acknowledged the limitation that our method may not be\nsuitable for organs without hierarchical organization, such as the liver, whose metabolic\nzonation 25 is continuously distributed along the porto-central axis 54.\nThe following revisions have been made in the manuscript:\n• Introduction: “Many tissues exhibit inherent hierarchical organization, wherein large\nstructures contain smaller, functionally specialized substructures. This pattern is evident\nin the layered architecture of the brain, the zonated morphology of lymphoid organs\nsuch as the spleen, and the compartmentalized ecosystem of tumors 25-27.”\n• Discussion: “As HRCHY-CytoCommunity explicitly models tissue organization as a\nhierarchy, it is particularly well-suited for organs with clearly stratified anatomical or\nfunctional layers. Its performance may be limited in tissues lacking strong hierarchical\norganization, such as liver, where metabolic zonation is continuously distributed\nalong the porto-central axis 54.”\nR3C2\n2. Evaluation on More Diverse Datasets: The manuscript focuses primarily on datasets with\nwell-separated cell types and pre-annotated structures. To evaluate the robustness of HRCHY-\nCytoCommunity, we recommend including sequencing-based spatial transcriptomics\ndatasets—especially those from the NeST study or similar benchmarks—to compare with\nmethods like NeST, which rely on transcriptome-wide data.\nWe thank the reviewer for this important suggestion. In the revised manuscript, we have\nexpanded the evaluation of HRCHY-CytoCommunity to include four representative\nsequencing-based spatial transcriptomics datasets covering different resolutions and\ntechnologies, including two datasets (the mouse hippocampus Slide-seq V2 dataset and the\nhuman breast cancer Visium V1 dataset) originally used in the NeST study.\n1. Mouse hippocampus Slide-seq V2 dataset (10 µm-spot resolution) HRCHY-\nCytoCommunity recovered biologically meaningful hierarchical structures spanning from\ncoarse-grained tissue compartments (TCs) including the cortical region, hippocampus, and\nbrain stem (Response Fig. 18b, middle panel) to fine-grained cellular neighborhoods (CNs)\nsuch as CA1-CA3 regions, dentate gyrus, corpus callosum, and distinct cortical\nlayers(Response Fig. 18c, middle panel), aligned well with the Allen Mouse Brain Atlas.\nIn contrast, NeST produced small and fragmented substructures rather than large-scale\ncoherent TCs at the coarse-grained level. Moreover, it left a considerable fraction of spots\nunassigned (labeled as “unidentified”) for both TC and CN levels, limiting the\ninterpretability of the tissue structure (Response Figs. 18b and c, right panels).\n2. Human breast cancer Visium V1 dataset (55 µm-spot resolution) HRCHY-\nCytoCommunity accurately distinguished tumor and non-tumor compartments (Response\nFig. 19b, middle panel) and recovered CNs closely matching manual annotations (Response\nFig. 19c, middle panel). While NeST correctly identified several CNs consistent with\nmanual annotations (e.g., DCIS/LCIS_1, DCIS/LCIS_3, and IDC_4), it incorrectly\npartitioned discontinuous tumor regions into multiple fragmented TCs at the coarse-grained\nlevel (Response Fig. 19b; right panel). Similar to the Slide-seq V2 case, NeST left the most\nspots located in Healthy_1 CN unassigned, reducing interpretability and completeness of\nthe tissue structure (Response Fig. 19c; right panel).\n3. Human colorectal cancer (CRC) Visium HD dataset (8 µm-bin resolution) HRCHY-\nCytoCommunity accurately distinguished tumor and normal compartments across all\nsamples, with moderate inclusion of intestinal epithelial-enriched regions within the tumor\ncompartment in some samples (Response Fig. 20b, left column). Conversely, NeST\nincorrectly assigned most bins to the tumor compartment (Response Fig. 20b, right column).\nAt the fine-grained level, HRCHY-CytoCommunity identified spatially coherent CNs\nacross samples (Response Fig. 20c, left column), while NeST produced fragmented and less\ninterpretable CNs (Response Fig. 20c, right column, cyan circles). Quantitatively, HRCHY-\nCytoCommunity achieved significantly higher Macro-F1 scores than NeST for coarse-\ngrained TC identification across all three samples (one-sided paired t-test P-values < 0.01;\nResponse Fig. 20d, right panel) and substantially higher average AMI values than NeST\n(Response Fig. 20d, left panel).\n4. Mouse intracerebral hemorrhage (ICH) Stereo-seq dataset (single-cell resolution)\nHRCHY-CytoCommunity effectively detected biologically meaningful fine-grained CNs,\nsuch as olfactory areas (OLF), meninges, and pallidum (PAL) (Response Fig. 21b, left\ncolumn), whereas NeST failed to capture these CNs in most samples (Response Fig. 21b,\nright column). At the coarse-grained level, HRCHY-CytoCommunity consistently\nidentified three biologically meaningful compartments: TC-0 (brain stem), TC-1 (cortex),\nand TC-2 (striatum), while NeST failed to distinguish these TCs. (Response Fig. 21c)\nQuantitatively, HRCHY-CytoCommunity achieved significantly higher AMI and Macro-\nF1 scores than NeST across all nine samples in fine-grained CN identification (one-sided\npaired t-test P-values < 1E-7; Response Fig. 21d).\nTaken together, these results demonstrate that HRCHY-CytoCommunity generalizes well to\nsequencing-based datasets and consistently recovers biologically interpretable hierarchical\ntissue structures across technologies and resolutions. Compared with NeST, which often\nproduces fragmented substructures and leaves many spots unassigned, HRCHY-\nCytoCommunity effectively captures both large- and small-scale spatially coherent tissue\nstructures and achieves complete cellular coverage.\nFull details of the extended performance evaluation on sequencing-based datasets are provided\nin the Results section under the titles: “Performance evaluation on coarse-grained TC\nidentification”, “Performance evaluation on fine-grained CN identification”, and “Cross-\nplatform generalization of HRCHY-CytoCommunity to low-resolution spatial\ntranscriptomics”, as well as in Supplementary Notes 1 and 2.\nResponse Fig. 18 (extracted from Fig. 4 and Supplementary Fig. 6). Performance evaluation of\nHRCHY-CytoCommunity on the mouse hippocampus Slide-seq V2 dataset. (a) Cluster stability measure\nused by HRCHY-CytoCommunity to determine the optimal numbers of coarse-grained TCs and fine-grained\nCNs. The Fowlkes-Mallows Index-based cluster stability values under each parameter setting are connected\nby solid lines. The red circle marks the selected stable resolution of tissue structures. (b) Coarse-grained TC\nreference from the Allen Mouse Brain Atlas 2 (left), and TCs identified by HRCHY-CytoCommunity (middle)\nand NeST (right) in the mouse hippocampus. (c) Fine-grained CN reference from the Allen Mouse Brain\nAtlas (left), and CNs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the mouse\nhippocampus. V3, third ventricle; MH, medial habenula; LH, lateral habenula; LP, lateral posterior nucleus\nof the thalamus; LD, lateral dorsal nucleus of the thalamus. “Unidentified” indicates cells not assigned to any\nTC or CN. “Unmatched” denotes TCs or CNs that could not be aligned with the Allen Mouse Brain Atlas .\nResponse Fig. 19 (extracted from Fig. 4 and Supplementary Fig. 7). Performance evaluation of\nHRCHY-CytoCommunity on the human breast cancer Visium V1 dataset. (a) Cluster stability measure\nused by HRCHY-CytoCommunity to determine the optimal numbers of coarse-grained TCs and fine-grained\nCNs. The Fowlkes-Mallows Index-based cluster stability values under each parameter setting are connected\nby solid lines. The red circle marks the selected stable resolution of tissue structures. (b) Manual annotation\nof coarse-grained TCs from Xu et al.3 (left), and TCs identified by HRCHY-CytoCommunity (middle) and\nNeST (right) in the human breast cancer. (c) Manual annotation of fine-grained CNs from Xu et al. (left),\nand CNs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the human breast cancer.\n“Unidentified” indicates cells not assigned to any TC or CN. Red arrowheads indicate the IDC_3 region\nmanually annotated by Xu et al., which was reclassified by both HRCHY-CytoCommunity and NeST as a\nnon-tumor CN.\nResponse Fig. 20 (extracted from Fig. 2 and Supplementary Fig. 3). Performance comparison of\nHRCHY-CytoCommunity and NeST on the human CRC Visium HD dataset. (a) Spatial maps of human\nCRC samples (8 μm bins) generated by Visium HD technology. Bins are colored by cell types (left) and\nmanually annotated coarse-grained structures (right). (b) Coarse-grained TCs identified by HRCHY-\nCytoCommunity (left) and NeST (right). (c) Fine-grained CNs identified by HRCHY-CytoCommunity (left)\nand NeST (right). (d) AMI and Macro-F1 scores calculated using manual compartment annotations. Each\npoint corresponds to the performance on an individual sample, with horizontal bars indicating the mean\nperformance across all samples. Points from the same sample are connected by grey dashed lines. P-values\nwere calculated using one-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.\n“Unidentified” indicates cells not assigned to any TC or CN. “Unmatched” denotes TCs that could not be\naligned with manual annotations.\nResponse Fig. 21 (extracted from Supplementary Figs. 5 and 6). Performance comparison of HRCHY-\nCytoCommunity and NeST on the mouse ICH Stereo-seq dataset. (a) Spatial maps showing the\ndistribution of cell types (left) and molecularly defined regions (right). (b) Fine-grained CNs identified by\nHRCHY-CytoCommunity (left) and NeST (right). (c) Coarse-grained TCs identified by HRCHY-\nCytoCommunity (left) and NeST (right). (d) AMI and Macro-F1 scores calculated using molecularly defined\nregions. Each point corresponds to the performance on an individual sample, with horizontal bars indicating\nthe mean performance across all samples. Points from the same samples are connected by grey dashed lines.\nP-values were calculated using a one-sided paired t-test. *, P-value < 0.05; **, P-value < 0.01; ***, P-value\n< 0.001. “Unidentified” indicates cells not assigned to any TC or CN. “Unmatched” denotes TCs that could\nnot be aligned with manual annotations. L2/3, cortical layer 2 and cortical layer 3; L4, cortical layer 4; L5,\ncortical layer 5; L6, cortical layer 6; Lesion, Lesion region; injured_STR, injured striatum region;\ninjured_cortex, injured cortex region; injured_PAL, injured pallidum region; DG, dentate gyrus; FT, fiber\ntracts; HY, hypothalamus; OLF, olfactory areas; MH, medial habenula; CA, Ammon’s horn; STRd, striatum\ndorsal region; STRv, striatum ventral region; TH, thalamus; VS, ventricular systems; isl, islm, islands of\nCalleja and major island of Calleja; PAL, pallidum; RT, reticular nucleus of the thalamus; CAsp, field\nAmmon’s horn, pyramidal layer; LSX, lateral septal complex; sAMY, striatum-like amygdala nuclei; CTXsp,\ncortical subplate.\nR3C3\n3. Limitations of Using Cell-Type Annotation Only: The exclusive reliance on cell-type labels\nposes a major limitation. High-quality annotations may not always be available, and annotation\nresolution varies by dataset. Furthermore, critical spatial structures may be defined not by cell\ntypes but by cell-cell interactions, gradients, or activity states. These limitations should be\nacknowledged, and potential solutions or future directions discussed. For instance, the author\ncould compare the results of two different levels of annotation (finer vs coarse).\nWe thank the reviewer for raising this valuable point. Indeed, the current implementation of\nHRCHY-CytoCommunity relies exclusively on cell-type annotations as input features, which\nraises concerns about its dependence on annotation quality and resolution that may vary across\ndatasets. This design choice was motivated by the broad availability of cell-type annotations\nacross diverse spatial omics datasets, which facilitates the generalizability of HRCHY-\nCytoCommunity across datasets. Moreover, it enhances the interpretability of the identified\nhierarchical tissue structures.\nTo address this concern, we systematically evaluated the robustness of HRCHY-\nCytoCommunity to both cell-type label inaccuracies and annotation resolution variations\nusing two datasets with complex cell-type compositions (i.e., the mouse hypothalamic preoptic\nregion MERFISH dataset and the human TNBC MIBI-TOF dataset), as detailed below:\n1. Robustness to cell-type label inaccuracies: We used the hierarchical tissue structures\nidentified from unperturbed cell-type annotations as the reference. Cell-type labels were\nthen randomly perturbed by 10%–50%, and the accuracy and consistency of the resulting\nhierarchical structures were compared with the reference. The Fowlkes-Mallows Index\n(FMI) was used to quantify the consistency between perturbed and reference structures\nacross different noise levels.\na) In the MERFISH dataset, the identified hierarchical structures remained highly\nconsistent until the label error rate exceeded 40% (Response Fig. 22a; average FMI >\n0.5 when label error rate = 40%). The accuracy of the fine-grained CNs (measured\nby AMI and Macro-F1 scores) showed only a moderate decline, with performance\nlargely preserved under reasonable noise levels (Response Figs. 22b and c)\nb) In the MIBI-TOF dataset, hierarchical structures remained highly consistent across\nperturbation rates (Response Fig. 23a; average FMI > 0.65 when label error rate =\n50%), The accuracy of the identified coarse-grained TCs (measured by the proportion\nof neoplastic and immune cells correctly assigned to their corresponding TCs) was\nlargely preserved under reasonable levels of annotation noise, with only a moderate\ndecrease as the label error rate increased (Response Fig. 23b).\n2. Robustness to cell-type annotation resolutions: We further assessed robustness to\nannotation resolution by comparing resulting structures obtained under different levels of\ncell-type granularity (i.e., 15 vs. 9 cell types in MERFISH, 17 vs. 11 cell types in TNBC\nMIBI-TOF). The FMI was used to quantify the similarity between identified structures.\nSpecifically,\na) HRCHY-CytoCommunity consistently recovered stable hierarchical tissue structures\nacross annotation resolutions, achieving average FMI scores of 0.673 in the\nMERFISH dataset (Response Figs. 22d-f) and 0.590 in the MIBI-TOF dataset across\nsamples (Response Figs. 23c-e).\nIn summary, these analyses demonstrate that HRCHY-CytoCommunity is strongly robust\nto both annotation noise and annotation resolution variations, alleviating concerns regarding\nstrong dependence on high-quality annotations and specific annotation resolution.\nFurthermore, we explicitly acknowledge the limitation that this design may not capture tissue\nstructures defined by continuous gradients, transient states, or functional interactions that are\nnot well represented by discrete cell-type labels. We have added this limitation to the\nDiscussion section and outlined future directions.\nThe corresponding discussion of the limitation of HRCHY-CytoCommunity is shown below.\n“Moreover, while the use of discrete cell-type annotations enhances the method’s\ngeneralizability across platforms, it may not fully capture structures defined by continuous\nmolecular gradients, transient states, or cell-cell communication signals. Future versions could\nincorporate multi-view learning to integrate such continuous features. In the meantime,\nmethods like NicheCompass 20 (focused on cellular communication-driven structures) and\nONTraC 55 (designed for spatially continuous niche trajectories) can complement HRCHY-\nCytoCommunity to provide a more comprehensive perspective on tissue organization.”\nFull details of the robustness assessment are provided in Supplementary Note 5 and\nSupplementary Figs. 19 and 20.\nResponse Fig. 22 (Supplementary Fig. 19). Robustness assessment of HRCHY-CytoCommunity to\nlabel inaccuracies and annotation resolution variations using the mouse hypothalamic preoptic region\nMERFISH dataset. (a) Hierarchical tissue structure similarity under increasing cell-type label perturbation\nrates (10%-50%), quantified by Fowlkes-Mallows Index (FMI)-based robustness score. Higher values\nindicate greater structural consistency with the unperturbed reference. (b-c) Fine-grained CN identification\naccuracy under label perturbations, measured by Adjusted Mutual Information (AMI) (b) and Macro-F1\nscores (c). Each point represents performance on one sample, with horizontal bars represent the mean\nperformance across n = 5 samples. Performances (points) on the same sample are connected by grey dashed\nlines. Mean values under each label error rate are connected by blue solid lines. Shaded areas denote standard\ndeviation. (d) Single-cell spatial maps showing the distribution of 15 cell types (left), coarse-grained TCs\n(middle) and fine-grained CNs (right) identified by HRCHY-CytoCommunity based on this high-resolution\nannotation. (e) Single-cell spatial maps showing the distribution of 9 cell types (left), coarse-grained TCs\n(middle) and fine-grained CNs (right) identified by HRCHY-CytoCommunity based on this low-resolution\nannotation. (f) Hierarchical tissue structure similarity between different cell-type annotation resolutions,\nquantified by FMI-based robustness score. The red dashed line indicates the mean similarity (0.673) relative\nto the 15-cell-type annotation-based reference.\nResponse Fig. 23 (Supplementary Fig. 20). Robustness assessment of HRCHY-CytoCommunity to\nlabel inaccuracies and annotation resolution variations using the human TNBC MIBI-TOF dataset.\n(a) Hierarchical tissue structure similarity under increasing cell-type label perturbation rates (10%-50%),\nquantified by Fowlkes-Mallows Index (FMI)-based robustness score. Higher values indicate greater\nstructural consistency with the unperturbed reference. (b) Coarse-grained TC identification accuracy under\nlabel perturbations, measured by proportion of neoplastic and immune cells correctly assigned to their\ncorresponding TCs. Each point represents performance on one compartmentalized tumor sample, with\nhorizontal bars represent the mean performance across n = 15 samples. Performances (points) on the same\nsample are connected by grey dashed lines. Mean values under each label error rate are connected by blue\nsolid lines. Shaded areas denote standard deviation. (c) Single-cell spatial maps showing the distribution of\n17 cell types (left), coarse-grained TCs (middle) and fine-grained CNs (right) identified by HRCHY-\nCytoCommunity based on this high-resolution annotation. (d) Single-cell spatial maps showing the\ndistribution of 11 cell types (left), coarse-grained TCs (middle) and fine-grained CNs (right) identified by\nHRCHY-CytoCommunity based on this low-resolution annotation. (e) Hierarchical tissue structure\nsimilarity between different cell-type annotation resolutions, quantified by FMI-based robustness score.\nThe red dashed line indicates the mean similarity (0.590) relative to the 17-cell-type annotation-based\nreference.\nR3C4\n4. Comparison with Additional Methods: The current comparisons omit several state-of-the-art\nspatial domain identification tools, including SpaGCN, NicheCompass, CellCharter, Gaston-\nMix among others. While these do not explicitly identify hierarchical domains, many can be\nrun iteratively or parameter-tuned to extract structures at different scales. Including such\ncomparisons would provide a fairer performance benchmark.\nWe appreciate the reviewer for the constructive suggestion. Following the suggestion from\nReviewers 2 and 3, we have benchmarked HRCHY-CytoCommunity against six state-of-the-\nart non-hierarchical spatial domain detection methods, including SpaGCN, NicheCompass,\nCellCharter, SpaSEG, Giotto Suite, and GraphST, in the revised manuscript. Since these\nmethods do not explicitly identify hierarchical tissue structures, we applied them separately at\ncoarse-grained TC and fine-grained CN resolutions. Performances were evaluated using\nAdjusted Mutual Information (AMI) and Macro-F1 scores across two datasets with ground-\ntruth TCs and two datasets with ground-truth CNs.\nOur results demonstrated that HRCHY-CytoCommunity consistently achieved superior or\ncompetitive performance than non-hierarchical methods. Specifically:\n1. For coarse-grained TC identification:\na) On the mouse spleen CODEX dataset (imaging-based), HRCHY-CytoCommunity\nwas the only method that consistently resolved the biologically meaningful lymphoid\ncompartment containing B-zone, periarteriolar lymphoid sheath, and marginal zone\nacross all samples (Response Fig. 24b, left column), whereas non-hierarchical methods\nlargely failed to capture this hierarchical structure (GraphST succeeded only in one\nsample). Quantitatively, HRCHY-CytoCommunity achieved significantly higher AMI\nand Macro-F1 scores than most non-hierarchical methods (one-sided paired t-test P-\nvalues < 0.05), while exhibiting performance comparable to GraphST (Response Fig.\n24d).\nb) On the human colorectal cancer (CRC) Visium HD dataset (sequencing-based),\nHRCHY-CytoCommunity was the only method that accurately distinguished tumor\nand normal compartments in the P2CRC sample, while all non-hierarchical methods\nfailed to do so. Similarly, in the remaining samples, these two compartments were also\nclearly delineated, though with moderate inclusion of intestinal epithelial-enriched\nregions within the tumor compartment (Response Fig. 25b, left column).\nQuantitatively, HRCHY-CytoCommunity achieved the highest average AMI and\nMacro-F1 scores among all representative non-hierarchical methods (Response Fig.\n25d).\n2. For fine-grained CN identification:\na) On the mouse hypothalamic preoptic region MERFISH dataset (imaging-based),\nHRCHY-CytoCommunity successfully recovered CNs consistent with manually\ndelineated nuclei with clear boundaries (Response Fig. 26b, left column), whereas\nother methods identified intermixed CNs and failed to capture smaller regions like PS\nand PaAP in most samples. Quantitatively, HRCHY-CytoCommunity significantly\noutperformed all non-hierarchical methods in both AMI and Macro-F1 metrics (one-\nsided paired t-test, p < 0.05), with the exception of CellCharter on the AMI score and\nGraphST on the Macro-F1 score, where the differences were not statistically\nsignificant (Response Fig. 26d)\nb) On the mouse intracerebral hemorrhage (ICH) Stereo-seq dataset (sequencing-\nbased), HRCHY-CytoCommunity achieved high accuracy in identifying fine-grained\nCNs, particularly excelling in detecting small-scale CNs such as olfactory areas (OLF),\nmeninges, and pallidum (PAL) (Response Fig. 27b, left column). Quantitatively,\nHRCHY-CytoCommunity achieved significantly higher AMI and Macro-F1 scores\nthan most non-hierarchical methods (one-sided paired t-test, P-values < 0.001), with\nperformance comparable to NicheCompass (Response Fig. 27d)\nTaken together, these results indicate that HRCHY-CytoCommunity is not only quantitatively\ncompetitive but also qualitatively superior, as it identifies biologically meaningful structures\nmissed by non-hierarchical methods. This highlights that the explicit hierarchical\nframework yields more accurate, stable, and interpretable multi-scale tissue organization\nacross diverse technologies.\nFull implementation details of these non-hierarchical spatial domain detection methods are\nprovided in the Methods section under the title “Running of published methods”. Full details\nof the performance evaluation are provided in the Results section under the title:\n“Performance evaluation on coarse-grained TC identification” and “Performance\nevaluation on fine-grained CN identification”, as well as in Supplementary Notes 1 and 2.\nResponse Fig. 24 (extracted from Fig. 2) Performance evaluation of HRCHY-CytoCommunity on the\nmouse spleen CODEX dataset. (a) Spatial maps of healthy mouse spleen samples generated by CODEX\ntechnology. Cells are colored by cell types (left) and manually annotated compartments (right), which serve\nas the ground truth for coarse-grained TCs. (b) Coarse-grained TCs identified by HRCHY-CytoCommunity\n(left) and NeST (right). (c) Coarse-grained TCs identified by non-hierarchical methods, including CellCharter,\nNicheCompass, GraphST, and SpaGCN. (d) Adjusted Mutual Information (AMI) and Macro-F1 scores\ncalculated using manual compartment annotations. Each point corresponds to the performance on an\nindividual sample, with horizontal bars indicating the mean performance across all samples. Points from the\nsame sample are connected by grey dashed lines. P-values were calculated using one-sided paired t-tests. *,\nP-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. “Unidentified” indicates cells not assigned to any\nTC. “Unmatched” denotes TCs that could not be aligned with manual annotations.\nResponse Fig. 25 (extracted from Fig. 2). Performance evaluation of HRCHY-CytoCommunity on the\nhuman CRC Visium HD dataset. (a) Spatial maps of human CRC samples (8 μm bins) generated by Visium\nHD technology. Bins are colored by cell types (left) and manually annotated coarse-grained structures (right).\n(b) Coarse-grained TCs identified by HRCHY-CytoCommunity (left) and NeST (right). (c) Coarse-grained\nTCs identified by non-hierarchical methods, including CellCharter, NicheCompass, and SpaSEG. (d) AMI\nand Macro-F1 scores calculated using manual compartment annotations. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same sample are connected by grey dashed lines. P-values were calculated using\none-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. “Unidentified” indicates\ncells not assigned to any TC. “Unmatched” denotes TCs that could not be aligned with manual annotations.\nResponse Fig. 26 (Fig. 3). Performance evaluation of HRCHY-CytoCommunity on the mouse\nhypothalamic preoptic region MERFISH dataset. (a) Spatial maps of mouse hypothalamic preoptic\nregion samples generated by MERFISH technology. Cells are colored by cell types (left) and manually\nannotated hypothalamic nuclei (right), which were generated based on the outlines of hypothalamic nuclei\n(middle) and serve as the ground truth for fine-grained CNs. (b) Fine-grained CNs identified by HRCHY-\nCytoCommunity (left) and NeST (right). (c) Fine-grained CNs identified by non-hierarchical methods,\nincluding CellCharter, NicheCompass, GraphST, SpaGCN, SpaSEG, and Giotto Suite. (d) Adjusted Mutual\nInformation (AMI) and Macro-F1 scores calculated using manually annotated hypothalamic nuclei. Each\npoint corresponds to the performance on an individual sample, with horizontal bars indicating the mean\nperformance across all samples. Points from the same sample are connected by grey dashed lines. P-values\nwere calculated using one-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.\n“Unidentified” indicates cells not assigned to any CN.\nResponse Fig. 27 (Supplementary Fig. 5). Performance evaluation of HRCHY-CytoCommunity on the\nmouse ICH Stereo-seq dataset. (a) Spatial maps showing the distribution of cell types (left) and molecularly\ndefined regions (right). (b) Fine-grained CNs identified by HRCHY-CytoCommunity (left) and NeST (right).\n(c) Fine-grained CNs identified by non-hierarchical tissue structure identification methods, including\nCellCharter, NicheCompass, GraphST, and SpaGCN. “Unidentified” indicates cells not assigned to any CN.\n(d) AMI and Macro-F1 scores calculated using molecularly defined regions. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same samples are connected by grey dashed lines. P-values were calculated using\na one-sided paired t-test. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. L2/3, cortical layer 2\nand cortical layer 3; L4, cortical layer 4; L5, cortical layer 5; L6, cortical layer 6; Lesion, Lesion region;\ninjured_STR, injured striatum region; injured_cortex, injured cortex region; injured_PAL, injured pallidum\nregion; DG, dentate gyrus; FT, fiber tracts; HY, hypothalamus; OLF, olfactory areas; MH, medial habenula;\nCA, Ammon’s horn; STRd, striatum dorsal region; STRv, striatum ventral region; TH, thalamus; VS,\nventricular systems; isl, islm, islands of Calleja and major island of Calleja; PAL, pallidum; RT, reticular\nnucleus of the thalamus; CAsp, field Ammon’s horn, pyramidal layer; LSX, lateral septal complex; sAMY,\nstriatum-like amygdala nuclei; CTXsp, cortical subplate.\nR3C5\n5. Loss Threshold Parameter (β): The parameter β is introduced as a threshold for selecting\nstructures based on graph loss but is not reflected in the mathematical formulations (Equations\n7–9). The manuscript should be updated to incorporate β directly into the equations, with a\nformal explanation of its role and its effect on structure selection.\nWe thank the reviewer for this valuable feedback. In the original manuscript, the loss threshold\nparameter β was introduced as a heuristic strategy to mitigate the cluster collapse problem (i.e.,\nall cells being assigned to a single cluster). In the revised version, we have replaced this\nheuristic approach with a trainable entropy-based balance regularization term to solve\nthis issue, which is now directly incorporated into the optimization objective. The formulation\nis defined as:\n𝑀𝑀 (2)\n1 𝐻𝐻𝑚𝑚\n𝐿𝐿𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 = 1− \n𝑀𝑀𝑚𝑚=1log𝑐𝑐2\nwhere is the number of perturbations, is the number of coarse-grained TCs. is the\n(2)\n𝑀𝑀 𝑐𝑐2 𝐻𝐻𝑚𝑚\nentropy of the cluster assignment distribution for the m-th perturbation. The is computed\n(2)\nas: 𝐻𝐻𝑚𝑚\n𝑏𝑏2\n(2) (2) (2)\n𝐻𝐻𝑚𝑚 = −𝑠𝑠𝑚𝑚,𝑘𝑘log𝑠𝑠𝑚𝑚,𝑘𝑘\n𝑘𝑘=1\nwhere represents the assignment probability of coarse-grained TC k in the m-th\n(2)\nperturba𝑠𝑠t𝑚𝑚io,n𝑘𝑘. This regularization term encourages balanced cluster assignments by maximizing\nthe entropy of the predicted distributions across clusters, thereby effectively preventing collapse\nof the coarse-grained TC assignments.\nAccordingly, the original β threshold is no longer used. Its role has been fully replaced by\nthe entropy-based regularization, which provides a more rigorous and differentiable solution\nto the same problem.\nFull details of this entropy-based balance regularization term are provided in the Methods\nsection under the title “Robust hierarchical tissue structure identification” (Equations 17-\n18).\nR3C6\n6. Details of Model Parameters: Please add a table summarizing the hyperparameters used for\neach dataset: number of GNN layers, MLP layers, values for α, β, and pruning thresholds.\nWithout this information, reproducibility and understanding of model performance across\ndatasets are limited.\nWe appreciate the reviewer’s constructive suggestion. In the revised manuscript, we have added\na comprehensive table (Supplementary Table 2) summarizing the hyperparameters used in all\ndatasets, and also described the key information below for clarity. Notably, in the revised\nversion, the original fixed α that balanced the unsupervised losses of the fine-grained CN and\nthe coarse-grained TC identification has been replaced by an adaptive α scheduling strategy\n(Methods), the previous orthogonality loss threshold β has been replaced by a trainable\nbalance regularization term (Methods), and the fixed edge-pruning threshold has been\nreplaced by an adaptive edge-pruning strategy (Methods).\nAccordingly, the following parameters are included in the table:\n• #TC and #CN, the numbers of TCs and CNs, respectively;\n• #K, the number of neighbors K in the KNN graph (cell-cell proximity graph)\nconstruction;\n• #GNN layers and #MLP layers, the numbers of the graph neural network and\nmultilayer perceptron modules, respectively;\n• #Epochs, the number of training epochs;\n• λ_consis and λ_balance, the weights of consistency and balance regularization terms,\nrespectively;\n• #Perturbations, the number of perturbations on the cell-cell proximity graph;\n• DropNode rate, the fraction of nodes whose features are masked during consistency\nregularization training;\n• Learning rate, the learning rate of the training process.\nDrop\n#GNN #MLP #Perturb Learnin\nDataset #TC #CN #K #Epochs Node\nlayers layers ations g rate\nrate\n𝝀𝝀𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 𝝀𝝀𝒃𝒃𝒃𝒃𝒃𝒃𝒃𝒃𝒄𝒄𝒄𝒄𝒃𝒃\nHuman CRC\nVisium HD 2 10 20\ndataset\nHuman Breast\nCancer Visium 2 9 20\ndataset\nMouse ICH\nStereo-seq 3 20 50\ndataset\nMouse spleen\nCODEX 2 7 50\ndataset\n2 2 1500 1 1 5 0.5 0.0001\nMouse\nhippocampus\n3 11 50\nSlide-seq V2\ndataset\nMouse\nhypothalamic\npreoptic region 2 12 100\nMERFISH\ndataset\nHuman TNBC\nMIBI-TOF 2 14 50\ndataset\nHuman Breast\nCancer IMC 2 7 20\ndataset\nSupplementary Table 2. Hyperparameter settings of HRCHY-CytoCommunity used in all datasets.\nThis table details the configuration of the HRCHY-CytoCommunity model for each dataset, including the\nnumber of tissue compartments (TCs) and cellular neighborhoods (CNs), the value of K in the K-nearest\nneighbor graph, the number of GNN and MLP layers, total training epochs, regularization weights (α,\nλ_consis, and λ_balance), the number of perturbations on the cell-cell proximity graph, DropNode rate, and\nlearning rate.\nR3C7\n7. Sensitivity Analysis: The study would benefit from conducting and reporting a sensitivity\nanalysis for α, β, and edge pruning on at least one image-based and one sequencing-based\ndataset. This is critical to demonstrate model stability and the impact of key hyperparameters\non performance.\nWe thank the reviewer for this valuable suggestion. In the revised manuscript, we have\nconducted comprehensive sensitivity analyses of key hyperparameters in HRCHY-\nCytoCommunity on both simulated and real datasets. Specifically, for coarse-grained TC\nidentification, we performed sensitivity analysis on the mouse spleen CODEX dataset\n(imaging-based) and the human CRC Visium HD dataset (sequencing-based). For fine-\ngrained CN identification, we performed the analysis on the mouse hypothalamic preoptic\nregion MERFISH dataset (imaging-based) and the mouse ICH Stereo-seq dataset\n(sequencing-based).\nAs mentioned in R3C6, in the revised manuscript, the original fixed α, β, and edge pruning\nthreshold have been replaced with an adaptive α scheduling strategy, a trainable entropy-based\nbalance regularization term, and an adaptive edge-pruning strategy, respectively. Therefore, our\nsensitivity analysis focused on five key hyperparameters: (1) the DropNode rate (fraction of\nnodes whose features are masked during consistency regularization training), (2) λ_consis\n(weight of the consistency regularization term), (3) λ_balance (weight of the balance\nregularization term), (4) the number of neighbors K in the KNN graph (cell-cell proximity\ngraph) construction, and (5) the number of perturbations used for generating perturbed cell-\ncell proximity graphs.\nThe results showed that:\n1. On the simulated dataset, HRCHY-CytoCommunity exhibited stable performance across a\nwide range of hyperparameter values, with only a slight decline in accuracy as the\nDropNode rate increased (Response Fig. 28).\n2. On real datasets,\na) For coarse-grained TC identification, performance was moderately sensitive to\nλ_consis and λ_balance, reflecting their roles in maintaining structural consistency and\npreventing cluster collapse. Excessively high DropNode rates degraded performance,\nlikely due to loss of informative features. In contrast, the adaptive edge-pruning\nstrategy ensured robustness to changes in the number of K neighbors. Increasing the\nnumber of perturbations did not substantially alter results (Response Fig. 29).\nb) For fine-grained CN identification, performance remained highly stable across all\nhyperparameter settings (Response Fig. 30).\nTaken together, these results demonstrate that HRCHY-CytoCommunity maintains stable and\naccurate performance without requiring extensive hyperparameter tuning.\nFull details of the sensitivity analysis are provided in Supplementary Note 4 and\nSupplementary Figs. 16-18.\nResponse Fig. 28 (Supplementary Fig. 16). Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity using the simulated dataset. Hierarchical tissue structure identification under different\nmodel configurations. Each configuration was run 10 times with different random seeds. Boxplot elements\nare defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range.\nPerformances are quantified by Adjusted Mutual Information (AMI) and Macro-F1 scores. Each point\ncorresponds to the performance of one run. This sensitivity analysis includes five hyperparameters: (1) the\nDropNode rate (fraction of nodes whose features are masked during consistency regularization training), (2)\nλ_consis (weight of the consistency regularization term), (3) λ_balance (weight of the balance regularization\nterm), (4) the number of neighbors K in the KNN graph (cell-cell proximity graph) construction, and (5) the\nnumber of perturbations used for generating perturbed cell-cell proximity graphs.\nResponse Fig. 29 (Supplementary Fig. 17). Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity on coarse-grained TC identification using real spatial omics datasets. Hierarchical\ntissue structure identification under different model configurations using (a) the mouse spleen CODEX\ndataset (imaging-based) and (b) the human colorectal cancer (CRC) Visium HD dataset (sequencing-based).\nEach run was run 10 times with different random seeds. Boxplot elements are defined as: center line, median;\nbox limits, upper and lower quartiles; whiskers, 1.5× interquartile range. Performances are quantified by\nAdjusted Mutual Information (AMI) and Macro-F1 scores. Each point corresponds to the performance of\none run. This sensitivity analysis includes five hyperparameters: (1) the DropNode rate (fraction of nodes\nwhose features are masked during consistency regularization training), (2) λ_consis (weight of the\nconsistency regularization term), (3) λ_balance (weight of the balance regularization term), (4) the number\nof neighbors K in the KNN graph (cell-cell proximity graph) construction, and (5) the number of perturbations\nused for generating perturbed cell-cell proximity graphs.\nResponse Fig. 30 (Supplementary Fig. 18). Sensitivity analysis of hyperparameters in HRCHY-\nCytoCommunity on fine-grained CN identification using real spatial omics datasets. Hierarchical tissue\nstructure identification under different model configurations using (a) the mouse hypothalamic preoptic\nregion MERFISH dataset (imaging-based) and (b) the mouse intracerebral hemorrhage (ICH) Stereo-seq\ndataset (sequencing-based). Each configuration was run 10 times with different random seeds. Boxplot\nelements are defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5×\ninterquartile range. Performances are quantified by Adjusted Mutual Information (AMI) and Macro-F1\nscores. Each point corresponds to the performance of one run. This sensitivity analysis includes five\nhyperparameters: (1) the DropNode rate (fraction of nodes whose features are masked during consistency\nregularization training), (2) λ_consis (weight of the consistency regularization term), (3) λ_balance (weight\nof the balance regularization term), (4) the number of neighbors K in the KNN graph (cell-cell proximity\ngraph) construction, and (5) the number of perturbations used for generating perturbed cell-cell proximity\ngraphs.\nR3C8\n8. Choice and Estimation of Number of Structures (c₁, c₂): The model requires the number of\nfine- and coarse-grained structures to be pre-defined. The manuscript should clarify whether\nthese are learned or fixed, and how users should choose them. If large values are set, how does\nthe model handle oversegmentation? A discussion or heuristic strategy would be\nhelpful.\nWe appreciate the reviewer’s constructive suggestions. In the revised manuscript, inspired by\nthe state-of-the-art spatial domain detection method CellCharter1, we have provided a heuristic\nstrategy to guide users in choosing the numbers of structures (c₁, c₂).\nIn this approach, users were first required to define a candidate search range for (c₁, c₂) and the\nnumber of repeated runs R. For each candidate pair (c₁, c₂), the model was trained R times\nindependently. The stability of the resulting hierarchical structures (i.e., CNs and TCs) was then\nquantified using the FMI computed across runs. Specifically, the FMI was computed separately\nfor the identified CNs and TCs across all pairs of runs for a given (c₁, c₂) configuration. This\nyielded a set of CN-level FMIs and a set of TC-level FMIs. For each level, the median FMI\nacross all run pairs was computed to summarize the cluster stability at that level. The final\nstability score for a (c , c ) configuration was defined as the average of the CN-level and TC-\n1 2\nlevel median FMIs. The optimal (c₁, c₂) pair was selected as the one that yields the highest\nstability score. Finally, among the R models trained under the optimal (c₁, c₂) configuration, the\nmodel achieving the lowest training loss was chosen to output the final hierarchical tissue\nstructures.\nThis procedure also addresses the reviewer’s concern regarding potential over-segmentation\nwhen large values of (c₁, c₂) are used:\n• When the number of clusters is too large, the clustering results tend to become unstable\nacross repeated runs, which is reflected by lower FMI stability.\n• Conversely, when the chosen (c₁, c₂) properly matches the underlying tissue\norganization, the model yields reproducible partitions with high FMI stability.\nIn summary, instead of requiring users to predefine (c₁, c₂), we introduce a data-driven\nheuristic strategy to determine appropriate values of (c₁, c₂) and prevent over-segmentation\nthrough stability analysis. This ensures that HRCHY-CytoCommunity can identify robust and\nstable hierarchical tissue structures.\nFull details are provided in the Methods section under the title “Determination of optimal\nnumbers of hierarchical tissue structures based on cluster stability”\n.\nR3C9\n9. Interpretability of Nested Structures: The manuscript suggests that assigning each fine-\ngrained structure to a single coarse-grained one improves interpretability. However, biological\nsystems often contain overlapping or ambiguous structures. For example, cells in a boundary\nregion may legitimately belong to multiple domains. Instead of enforcing strict assignments,\nthe model might benefit from allowing flexible memberships or highlighting ambiguity.\nWe thank the reviewer for this insightful comment. While HRCHY-CytoCommunity currently\nenforces strict hierarchical assignments for clear interpretability and consistency in downstream analyses, we\nfully agree that cells located in boundary regions may exhibit ambiguity and overlapping\nmemberships in hierarchical structures. To address this concern, we provided an\ninterface that allows users to visualize regions with potentially ambiguous hierarchical\nmemberships. Specifically, we introduced an uncertainty measure based on the entropy of the\nsoft assignment distribution for each cell. As shown in Response Fig. 31, cells with higher\nentropy scores are frequently located at CN boundaries or transition zones, reflecting potential\nambiguity in structural assignments. This interface enables users to assess assignment\nuncertainty and highlight areas where nested or overlapping structures may exist. In future work,\nwe plan to extend the framework by explicitly incorporating this uncertainty into probabilistic\nor soft-membership hierarchical models.\nFull details of the ambiguous structure analysis are provided in Supplementary Note 9.\nResponse Fig. 31 (Supplementary Fig. 24). Ambiguous structures in the mouse Slide-seq V2 dataset.\n(a) Fine-grained CNs identified by HRCHY-CytoCommunity. (b) Uncertainty score distribution, with values\nnear 0 indicating high confidence in tissue structure assignment and values near 1 reflecting high ambiguity.\nR3C10\n10. Survival Analysis Details: The survival analysis lacks details on the inclusion of clinical\ncovariates. Without covariate adjustment, survival differences may be confounded. If covariates\n(e.g., age, sex, tumor grade) were not included, either re-run the analysis with appropriate\nadjustments or justify why this analysis is still valid.\nWe thank the reviewer for bringing this critical detail up and fully agree that covariate\nadjustment is essential for survival analysis. Accordingly, in the revised manuscript, we have\nreplaced all previous results derived from log-rank tests with those obtained from\nmultivariable Cox proportional hazards models.\nIn this analysis, we included three well-established prognostic covariates for breast cancer (age,\ntumor grade, and tumor size). Sex was not included as a covariate because all patients in this\ndataset were female. The results were reported as hazard ratios (HRs) and P-values (e.g., HR =\n2.2, P-value = 0.031).\nOverall, our results suggest the potential of hierarchical tissue structures identified by HRCHY-\nCytoCommunity serving as informative features for patient stratification, and reveal subgroups\nwith independent prognostic significance.\nFull implementation details of the survival analysis are provided in the Methods section under\nthe title “Survival analysis”, and the corresponding results in the human breast cancer IMC\ndataset are provided in the Results section under the title “Hierarchical tissue structure-\nbased stratification of breast cancer patients”\n.\nR3C11\n11. Validation with Simulated Data: Consider including a synthetic dataset with known\nhierarchical ground truth to benchmark model accuracy and interpretability under controlled\nconditions.\nWe thank the reviewer for this constructive suggestion. Following the suggestion, we\nconstructed a synthetic dataset with a known hierarchical ground truth to benchmark\nHRCHY-CytoCommunity under controlled conditions. Specifically, inspired by the human\nTNBC MIBI-TOF dataset, we designed two TCs, one dominated by immune cells (cell type 1)\nand the other by tumor cells (cell type 8). Each TC contained multiple CNs with distinct\ncompositional patterns, yielding a ground-truth structure of two TCs and eight CNs. The cell-\ntype composition profiles for CNs in simulation data has been provided in Supplementary\nTable 4.\nOur results showed that HRCHY-CytoCommunity accurately recovered the ground-truth\nhierarchical structures, achieving AMI = 0.956 and Macro-F1 score = 0.994 for TCs, and\nAMI = 0.914 and Macro-F1 score = 0.959 for CNs (Response Fig. 32). We further leveraged\nthis dataset to perform ablation and sensitivity analyses. The ablation experiments confirmed\nthat the consistency and balance regularization terms are most effective when applied together\nand also validated the effectiveness of the adaptive edge-pruning strategy in coarse-grained TC\nidentification (Response Fig. 33). The sensitivity analysis demonstrated that model\nperformance remains stable across a wide range of hyperparameter choices (Response Fig. 28).\nTaken together, the experiments using the synthetic dataset verify that HRCHY-\nCytoCommunity can accurately identify hierarchical tissue structures under controlled\nconditions and further highlight the effectiveness and robustness of the proposed design.\nFull details of the construction of the simulated dataset are provided in Supplementary Note\n7, under the title “Simulated dataset”, and the results of HRCHY-CytoCommunity on this\ndataset are provided in Supplementary Fig. 23. Full details of the ablation and sensitivity\nanalyses using this dataset are provided in Supplementary Notes 3 and 4.\nResponse Fig. 32 (Supplementary Fig. 23). Performance evaluation of HRCHY-CytoCommunity using\nthe simulated dataset. (a) Spatial distribution of eight cell types in the simulated dataset. (b) Ground truth\nof the coarse-grained TCs. (c) Coarse-grained TCs identified by HRCHY-CytoCommunity (AMI = 0.956,\nMacro-F1 score = 0.994). (d) Ground truth of the fine-grained CNs. (e) Fine-grained CNs identified by\nHRCHY-CytoCommunity (AMI = 0.914, Macro-F1 score= 0.959).\nResponse Fig. 33 (Supplementary Fig. 13). Ablation study of key components in HRCHY-\nCytoCommunity using the simulated dataset. Hierarchical tissue structure identification under different\nmodel configurations. Each configuration was run 10 times with different random seeds. Boxplot elements\nare defined as: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range.\nPerformances are quantified by Adjusted Mutual Information (AMI) and Macro-F1 scores. Each point\ncorresponds to the performance of one run. Ablation studies were performed on (a) consistency and balance\nregularization terms, (b) the adaptive edge-pruning strategy, and (c) the adaptive α scheduling strategy.\nR3C12\n12. Multiple Sample Integration: The method currently runs HRCHY-CytoCommunity\nindependently per sample, followed by clustering enrichment profiles. This strategy might not\nfully capture shared spatial structure across samples. Consider discussing or exploring\nalternatives, such as joint embedding approaches or graph alignment techniques.\nWe thank the reviewer for this insightful comment. We agree that joint embedding or graph\nalignment strategies could better capture shared spatial structures across samples. To explore\nthis, we implemented a straightforward multi-sample integration by merging individual sample-\nspecific KNN graphs into a single large graph. In this integrated graph, each connected\ncomponent corresponded to a sample-specific KNN graph. We then applied HRCHY-\nCytoCommunity directly to this integrated multi-sample graph.\nThe model consistently identified symmetric CNs that were highly concordant with the\nmanually delineated nuclei (Response Fig. 34c). Quantitatively, the multi-sample strategy\nyielded a slightly lower AMI compared with the single-sample setting (two-sided paired t-test,\nP-value = 0.034; Response Fig. 34d, left panel) and a comparable Macro-F1 score (two-sided\npaired t-test, P-value = 0.55; Response Fig. 34d, right panel). These results suggested that this\nmulti-sample integration is feasible, although sample-specific heterogeneity may slightly\nreduce accuracy. However, given that HRCHY-CytoCommunity is trained on the entire graph,\nscaling this integration to large cohorts would increase memory and computational\nrequirements, thereby limiting its practicality for large-scale studies. Moreover, in the TNBC\nMIBI-TOF and the breast cancer IMC datasets, we have demonstrated that clustering cell-type\nenrichment profiles across samples can effectively capture both conserved and sample-specific\nTCs and CNs.\nTherefore, we retained the single-sample implementation and used the clustering-based\napproach to perform multi-sample integration in this study, which offers a favorable trade-off\nbetween robustness and scalability. In future work, we plan to explore graph sampling and\nalignment-based techniques to develop a scalable cross-sample extension of HRCHY-\nCytoCommunity. This discussion has been incorporated into the revised Discussion section.\n“Additionally, future extensions could incorporate graph sampling or alignment techniques to\nimprove cross-sample integration within an end-to-end analytical workflow.”\nResponse Fig. 34. Performance evaluation of HRCHY-CytoCommunity with multi-sample integration\non the mouse hypothalamic preoptic region MERFISH dataset. (a) Spatial maps of mouse hypothalamic\npreoptic region samples generated by MERFISH technology. Cells are colored by cell types (left) and\nmanually annotated hypothalamic nuclei (right), which were generated based on the outlines of hypothalamic\nnuclei (middle) and serve as the ground truth for fine-grained CNs. (b) Hierarchical tissue structure identified\nby the single-sample version of HRCHY-CytoCommunity, coarse-grained TCs (left) and fine-grained CNs\n(right). (c) Hierarchical tissue structure identified by the multi-sample version of HRCHY-CytoCommunity,\ncoarse-grained TCs (left) and fine-grained CNs (right). (d) Adjusted Mutual Information (AMI) and Macro-\nF1 scores calculated using manually annotated hypothalamic nuclei. Each point corresponds to the\nperformance on an individual sample, with horizontal bars indicating the mean performance across all\nsamples. Points from the same sample are connected by grey dashed lines. P-values were calculated using\none-sided paired t-tests. *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.\nMinor Comments\nR3C1\n1. Please add a legend to Figure 3a to clarify the visualized elements and improve readability.\nWe thank the reviewer for this helpful suggestion to improve the readability of our manuscript.\nIn the revised version, we have added a clear cell-type legend below Figure 3a to clarify the\nvisualized elements and enhance readability.\nR3C2\n2. The rationale for selecting the “3TC” annotation in the bregma mouse brain sample is not\nprovided. Please justify its use or clarify its biological relevance.\nWe thank the reviewer for pointing this out. In the original manuscript, we empirically set the\nnumber of coarse-grained TCs in the mouse hypothalamic preoptic region MERFISH dataset\nto three. In the revised version, we employed the cluster stability-based heuristic strategy to\ndetermine the optimal number of coarse-grained TCs and fine-grained CNs in this dataset (see\nR3C8). For performance benchmarking, we fixed the number of fine-grained CNs to 12\n(corresponding to the number of ground-truth CNs) and applied this strategy to determine the\noptimal number of coarse-grained TCs. The results showed that the optimal number of coarse-\ngrained TCs is two (Response Fig. 35). Accordingly, we have updated the identified hierarchical\nstructures under this configuration in the revised manuscript.\nResponse Fig. 35 (extracted from Supplementary Fig. 4). Cluster stability-based heuristic strategy to\nidentify the optimal number of coarse-grained TCs and fine-grained CNs on the mouse hypothalamic\npreoptic region MERFISH dataset. (a) Cluster stability measure used by HRCHY-CytoCommunity to\ndetermine the optimal numbers of coarse-grained TCs and fine-grained CNs. Each point represents the\nFowlkes-Mallows Index-based cluster stability for an individual sample. Horizontal bars indicate the mean\nacross all samples. Mean values under each parameter setting are connected by grey solid lines. The red circle\nmarks the selected stable resolution of tissue structures.\nR3C3\n3. For the brain example, it would be informative to overlay the results with reference\nannotations from the Allen Brain Atlas to validate spatial domains at both coarse and fine levels.\nWe thank the reviewer for this important suggestion. In the revised manuscript, we have\nvalidated the spatial domains identified by HRCHY-CytoCommunity against the reference\nanatomical annotations from the Allen Mouse Brain Atlas at both coarse- and fine-grained\nlevels, using three representative brain datasets: the mouse hypothalamic preoptic region\n(MERFISH), the mouse hippocampus (Slide-seq V2), and the mouse ICH (Stereo-seq).\n1. Coarse-grained validation:\nIn the Slide-seq V2 dataset, the identified coarse-grained TCs exhibited a strong\ncorrespondence with major anatomical regions in the Allen Mouse Brain Atlas,\nincluding the cortical region, hippocampus, and brain stem (Response Fig. 18c, middle\npanel). Similarly, in the Stereo-seq dataset, HRCHY-CytoCommunity consistently\nidentified three biologically meaningful compartments: TC-0 (brain stem), TC-1\n(cortex), and TC-2 (striatum), which closely matched the reference atlas (Response Fig.\n21c, left column).\n2. Fine-grained validation:\nAt the fine-grained level, HRCHY-CytoCommunity accurately resolved hippocampal\nsubstructures in the Slide-seq V2 dataset, including the CA1-CA3 regions, dentate gyrus,\ncorpus callosum, and cortical layers (Response Fig. 18d, middle panel). In the Stereo-\nseq dataset, it effectively delineated fine-grained CNs such as cortical layers, olfactory\nareas (OLF), meninges, and pallidum (PAL), all of which aligned well with the fine-\nlevel Allen Mouse Brain Atlas annotations (Response Fig. 21b, left column). Likewise,\nin the MERFISH dataset, the identified CNs corresponded to known hypothalamic\nnuclei (e.g., medial and lateral preoptic nuclei), displaying symmetric spatial\norganization consistent with the reference atlas (Response Fig. 26b, left column).\nIn summary, these results demonstrate that HRCHY-CytoCommunity not only achieves high\nquantitative accuracy but also recovers spatially and anatomically coherent domains that are\nhighly concordant with the Allen Mouse Brain Atlas across multiple scales, thereby supporting\nthe biological interpretability of the hierarchical tissue structures identified by HRCHY-\nCytoCommunity.\nR3C4\n4. In some examples, fine-grained domains appear more as transitional gradients than discrete\nclusters. This challenges the assumption of discrete, nested hierarchies and should be discussed.\nWe thank the reviewer for this insightful comment and fully agree that in certain biological\ncontexts, fine-grained domains may exhibit continuous or transitional gradients rather than\ndiscrete cluster boundaries. In the revised manuscript, we have acknowledged this as a\nlimitation of the current framework and clarified that HRCHY-CytoCommunity models tissue\norganization as a discrete hierarchy to enhance interpretability. Nevertheless, the use of soft\nprobabilistic assignments allows the model to partially capture gradual transitions between\nadjacent regions (see R3C9, Response Fig. 31).\nWe also note that future extensions could incorporate trajectory-based representations to better\ncharacterize continuous spatial gradients. For example, recent approaches such as ONTraC 26\nreconstruct spatially continuous niche trajectories and could complement discrete hierarchical\nmodeling by providing a continuous perspective on tissue organization.\nThe corresponding statement has been added to the revised Discussion section:\n“Moreover, while the use of discrete cell-type annotations enhances the method’s\ngeneralizability across platforms, it may not fully capture structures defined by continuous\nmolecular gradients, transient states, or cell-cell communication signals. Future versions could\nincorporate multi-view learning to integrate such continuous features. In the meantime,\nmethods like NicheCompass 20 (focused on cellular communication-driven structures) and\nONTraC 55 (designed for spatially continuous niche trajectories) can complement HRCHY-\nCytoCommunity to provide a more comprehensive perspective on tissue organization.”\nR3C5\n5. While this method builds on Cytocommunity, the distinction lies primarily in the added CNN\nfor hierarchical domain inference. Please clarify what conceptual or performance advancements\nare made beyond the original Cytocommunity approach.\nWe thank the reviewer for raising this important question. While HRCHY-CytoCommunity\nbuilds upon the conceptual foundation of CytoCommunity, it introduces substantial\nconceptual, methodological, and performance advancements that enable hierarchical\ndomain inference, which is not supported by the original CytoCommunity framework.\nSpecifically:\n1. Conceptual advancement: CytoCommunity performs flat (non-hierarchical) tissue\nstructure detection, identifying only a single-scale community structure. In contrast,\nHRCHY-CytoCommunity reformulates the problem as a hierarchical community\ndetection task, jointly optimizing fine-grained CNs and coarse-grained TCs within a\nunified, fully differentiable framework. This design enables the model to capture multi-\nscale spatial organization that more faithfully reflects true tissue architecture.\n2. Methodological innovation: HRCHY-CytoCommunity introduces a two-level\ndifferentiable pooling mechanism based on graph minimum-cut optimization, enabling\nlearnable aggregation of fine-level CNs into coarse-level TCs. In addition, two new\nregularization strategies are incorporated:\na) Consistency regularization, which improves stability and reproducibility\nacross stochastic perturbations.\nb) Entropy-based balance regularization, which prevents cluster collapse and\nencourages balanced hierarchical assignments.\nIn summary, these conceptual and methodological advancements make HRCHY-\nCytoCommunity a distinct framework that advances the original CytoCommunity from flat\nclustering to hierarchical tissue structure identification.\nR3C6\n6. The claim that NeST is not applicable to low-feature datasets needs further explanation.\nSince NeST uses gene expression, it may still be applicable depending on the resolution. If\nexcluded from comparisons, the manuscript should clearly justify why.\nWe thank the reviewer for this valuable comment and apologize for any confusion caused by\nour previous description of NeST’s applicability to low-feature datasets. In the revised\nmanuscript, we have clarified the description of NeST to more accurately reflect its scope and\nlimitations:\n“However, NeST relies on gene expression features, limiting its performance on data with\nsparse molecular measurements, such as spatial proteomics datasets”\nMoreover, in the revised manuscript, we have expanded the performance benchmarking of\nHRCHY-CytoCommunity against NeST to include both imaging-based and sequencing-based\ndatasets spanning a wide range of feature dimensionalities, as listed below:\n1. Mouse spleen CODEX dataset (single-cell resolution, 30 proteins)\n2. Human CRC Visium HD dataset (8 µm-bin resolution, 18085 genes)\n3. Mouse hypothalamic preoptic region MERFISH dataset (single-cell resolution, 155 genes)\n4. Mouse ICH Stereo-seq dataset (single-cell resolution, 23625 genes)\n5. Mouse hippocampus Slide-seq V2 dataset (10 µm-spot resolution, 17733 genes)\n6. Human breast cancer Visium V1 dataset (55 µm-spot resolution, 24923 genes)\nAcross all datasets, HRCHY-CytoCommunity consistently achieved substantially superior\nperformance compared to NeST, both quantitatively (higher AMI and Macro-F1 scores)\nand qualitatively (more spatially coherent and biologically interpretable hierarchical structures)\n(Response Figs. 18-21, 24, 26).\nFull details of performance benchmarking are provided in the Results section under titles:\n“Performance evaluation on coarse-grained TC identification”, “Performance evaluation\non fine-grained CN identification”, and “Cross-platform generalization of HRCHY-\nCytoCommunity to low-resolution spatial transcriptomics”, as well as in Supplementary\nNotes 1 and 2.\nR3C7\n7. The use of manual annotations in unorganized organs (e.g., lung, kidney) is potentially\nproblematic, as these organs lack strong spatial hierarchy. Please explain the rationale behind\nthese choices and whether such annotations are used as training or just evaluation.\nWe thank the reviewer for this valuable comment. We fully agree that investigating hierarchical\ntissue structures in organs with unorganized or weak spatial hierarchies (e.g., lung or kidney)\ncan be problematic. Accordingly, in this study, HRCHY-CytoCommunity was applied only\nto organs with well-established multi-scale hierarchical organization that have been\nconsistently reported in previous studies, including lymphoid organs 10,12 (e.g., spleen),\ntumor tissues 15,17,30 (e.g., breast cancer and colorectal cancer), and the brain (e.g., hypothalamic\npreoptic region and hippocampus).\nTo avoid misunderstanding, we have clarified the scope of our method in the revised\nmanuscript. The following revisions have been made in the revised Discussion section:\n“As HRCHY-CytoCommunity explicitly models tissue organization as a hierarchy, it is\nparticularly well-suited for organs with clearly stratified anatomical or functional layers. Its\nperformance may be limited in tissues lacking strong hierarchical organization, such as liver,\nwhere metabolic zonation is continuously distributed along the porto-central axis 54.”\nReviewer #3 (Remarks on code availability):\nThe repository provided is well detailed and includes all information needed to reproduce the\nanalysis presented and use the tool.\nWe thank the reviewer for the positive feedback and are glad to hear that the repository was\nfound to be well documented and reproducible. The updated source code and documentation of\nHRCHY-CytoCommunity are available at https://github.com/huBioinfo/HRCHY-\nCytoCommunity.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This\nis part of the Nature Communications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe thank the co-reviewer for their time and valuable contribution to the peer-review process.\nReferences:\n1. Varrone, M., Tavernari, D., Santamaria-Martínez, A., Walsh, L. A. & Ciriello, G.\nCellCharter reveals spatial cell niches associated with tissue remodeling and cell plasticity.\nNat Genet 56, 74–84 (2024).\n2. Wang, Q. et al. The Allen Mouse Brain Common Coordinate Framework: A 3D Reference\nAtlas. Cell 181, 936-953.e20 (2020).\n3. Xu, H. et al. Unsupervised spatially embedded deep representation of spatial\ntranscriptomics. Genome Medicine 16, 12 (2024).\n4. Feng, W. et al. Graph Random Neural Networks for Semi-Supervised Learning on Graphs.\n5. Bai, Y. et al. SpaSEG: unsupervised deep learning for multi-task analysis of spatially\nresolved transcriptomics. Genome Biol 26, 230 (2025).\n6. Chen, J. G. et al. Giotto Suite: a multiscale and technology-agnostic spatial multiomics\nanalysis ecosystem. Nat Methods 1–13 (2025) doi:10.1038/s41592-025-02817-w.\n7. Long, Y. et al. Spatially informed clustering, integration, and deconvolution of spatial\ntranscriptomics with GraphST. Nat Commun 14, 1155 (2023).\n8. Hu, J. et al. SpaGCN: Integrating gene expression, spatial location and histology to identify\nspatial domains and spatially variable genes by graph convolutional network. Nat Methods\n18, 1342–1351 (2021).\n9. Birk, S. et al. Quantitative characterization of cell niches in spatially resolved omics data.\nNat Genet https://doi.org/10.1038/s41588-025-02120-6 (2025) doi:10.1038/s41588-025-\n02120-6.\n10. Mebius, R. E. & Kraal, G. Structure and function of the spleen. Nat Rev Immunol 5, 606–\n616 (2005).\n11. Cerutti, A., Cols, M. & Puga, I. Marginal zone B cells: virtues of innate-like antibody-\nproducing lymphocytes. Nat Rev Immunol 13, 118–132 (2013).\n12. Goltsev, Y. et al. Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed\nImaging. Cell 174, 968-981.e15 (2018).\n13. Lake, B. B. et al. An atlas of healthy and injured cell states and niches in the human kidney.\nNature 619, 585–594 (2023).\n14. Fridman, W. H. et al. B cells and tertiary lymphoid structures as determinants of tumour\nimmune contexture and clinical outcome. Nat Rev Clin Oncol 19, 441–457 (2022).\n15. Keren, L. et al. A Structured Tumor-Immune Microenvironment in Triple Negative Breast\nCancer Revealed by Multiplexed Ion Beam Imaging. Cell 174, 1373-1387.e19 (2018).\n16. Ohara, K. et al. The evolution of metastatic upper tract urothelial carcinoma through\ngenomic-transcriptomic and single-cell protein markers analysis. Nat Commun 15, 2009\n(2024).\n17. Oliveira, M. F. D. et al. High-definition spatial transcriptomic profiling of immune cell\npopulations in colorectal cancer. Nat Genet https://doi.org/10.1038/s41588-025-02193-3\n(2025) doi:10.1038/s41588-025-02193-3.\n18. Arora, R. et al. Spatial transcriptomics reveals distinct and conserved tumor core and edge\narchitectures that predict survival and targeted therapy response. Nat Commun 14, 5029\n(2023).\n19. Yan, Y. et al. Multi-omic profiling highlights factors associated with resistance to immuno-\nchemotherapy in non-small-cell lung cancer. Nat Genet 57, 126–139 (2025).\n20. MacFawn, I. P. et al. The activity of tertiary lymphoid structures in high grade serous\novarian cancer is governed by site, stroma, and cellular interactions. Cancer Cell 42, 1864-\n1881.e5 (2024).\n21. Meylan, M. et al. Tertiary lymphoid structures generate and propagate anti-tumor antibody-\nproducing plasma cells in renal cell cancer. Immunity 55, 527-541.e5 (2022).\n22. Tsutsumida, H. et al. A micropapillary pattern is predictive of a poor prognosis in lung\nadenocarcinoma, and reduced surfactant apoprotein A expression in the micropapillary\npattern is an excellent indicator of a poor prognosis. Mod Pathol 20, 638–647 (2007).\n23. Sautès-Fridman, C., Petitprez, F., Calderaro, J. & Fridman, W. H. Tertiary lymphoid\nstructures in the era of cancer immunotherapy. Nat Rev Cancer 19, 307–325 (2019).\n24. Schumacher, T. N. & Thommen, D. S. Tertiary lymphoid structures in cancer. Science\nhttps://doi.org/10.1126/science.abf9419 (2022) doi:10.1126/science.abf9419.\n25. Kietzmann, T. Metabolic zonation of the liver: The oxygen gradient revisited. Redox Biol\n11, 622–630 (2017).\n26. Wang, W., Zheng, S., Shin, S. C., Chávez-Fuentes, J. C. & Yuan, G.-C. ONTraC\ncharacterizes spatially continuous variations of tissue microenvironment through niche\ntrajectory analysis. Genome Biology 26, 117 (2025).\n27. Bhate, S. S., Barlow, G. L., Schürch, C. M. & Nolan, G. P. Tissue schematics map the\nspecialization of immune tissue motifs and their appropriation by tumors. Cell Syst 13, 109-\n130.e6 (2022).\n28. Hickey, J. W. et al. Organization of the human intestine at single-cell resolution. Nature\n619, 572–584 (2023).\n29. Nettekoven, C. et al. A hierarchical atlas of the human cerebellum for functional precision\nmapping. Nat Commun 15, 8376 (2024).\n30. Jackson, H. W. et al. The single-cell pathology landscape of breast cancer. Nature 578, 615–\n620 (2020).\nResponse to reviewer comments\nReviewer #1 (Remarks to the Author):\nIdentifying mulit-scale tissue structures from ST data is crucial for downstream\nbiological discovery. This work provides a novel method to identify hierarchical\nstructures.\nThis submission has made substantial revisions and addressed all my proposed\nquestions. I recommend acceptance of this work.\nReviewer #2 (Remarks to the Author):\nI am satisfied with the authors’ response to my review comments. The revisions are\ndetailed and directly address all major concerns raised in the initial review.\nSpecificially, the authors have substantially strengthened the work by expanding the\nbenchmarking, adding comprehensive ablation and robustness analyses, and most\nimportantly replacing the post hoc majority-voting strategy with a principled, end-to-\nend regularization framework. These changes not only improve the technical rigor of\nthe method but also enhance its transparency and practical relevance. I have no further\nconcerns.\nReviewer #3 (Remarks to the Author):\nThe authors have substantially improved the manuscript in comparison to the first\nversion. We are satisfied with the answers given to our concerns and with the\nmodifications done in text and figures. We would like to highlight the level of detail\nprovided in the answer’s to reviewer’s document, which better contextualize the\nbenefits and limitations of this method in comparison to previous ones. We\nrecommend the acceptance of the study.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed\nreports. This is part of the Nature Communications initiative to facilitate training in\npeer review and to provide appropriate recognition for Early Career Researchers who\nco-review manuscripts.\nWe sincerely thank all reviewers for the positive and encouraging comments\nregarding our revised manuscript. We are also grateful for their careful evaluation and\nconstructive feedback throughout the entire review process, which have significantly\nstrengthened our work."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Identifying mulit-scale tissue structures from ST data is crucial for downstream biological discovery. This work provides a\nnovel method to identify hierarchical structures.\nThis submission has made substantial revisions and addressed all my proposed questions. I recommend acceptance of this\nwork.", "reviewer_2": "I am satisfied with the authors’ response to my review comments. The revisions are detailed and directly address all major\nconcerns raised in the initial review. Specificially, the authors have substantially strengthened the work by expanding the\nbenchmarking, adding comprehensive ablation and robustness analyses, and most importantly replacing the post hoc\nmajority-voting strategy with a principled, end-to-end regularization framework. These changes not only improve the\ntechnical rigor of the method but also enhance its transparency and practical relevance. I have no further concerns", "reviewer_3": "The authors have substantially improved the manuscript in comparison to the first version. We are satisfied with the answers\ngiven to our concerns and with the modifications done in text and figures. We would like to highlight the level of detail\nprovided in the answer's to reviewer's document, which better contextualize the benefits and limitations of this method in\ncomparison to previous ones. We recommend the acceptance of the study.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts."}, "author_response": "We sincerely thank all reviewers for the positive and encouraging comments\nregarding our revised manuscript. We are also grateful for their careful evaluation and\nconstructive feedback throughout the entire review process, which have significantly\nstrengthened our work."}]

Decision letter