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

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Nat Commun . 2026 Feb 28;17:3312. doi: 10.1038/s41467-026-70069-z Search in PMC Search in PubMed View in NLM Catalog Add to search HRCHY-CytoCommunity identifies hierarchical tissue organization in cell-type spatial maps Runzhi Xie Runzhi Xie 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Runzhi Xie 1, # , Zekun Wang Zekun Wang 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Zekun Wang 1, # , Jianrui Liu Jianrui Liu 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Jianrui Liu 1 , Han Xu Han Xu 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Han Xu 1 , Yafei Xu Yafei Xu 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Yafei Xu 1 , Jiadong Lin Jiadong Lin 2 School of Automation Science and Engineering, Faculty of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi China 3 MOE Key Lab for Intelligent Networks & Networks Security, Faculty of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi China Find articles by Jiadong Lin 2, 3, ✉ , Yuxuan Hu Yuxuan Hu 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Yuxuan Hu 1, ✉ , Lin Gao Lin Gao 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China Find articles by Lin Gao 1, ✉ Author information Article notes Copyright and License information 1 School of Computer Science and Technology, Xidian University, Xi’an, Shaanxi China 2 School of Automation Science and Engineering, Faculty of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi China 3 MOE Key Lab for Intelligent Networks & Networks Security, Faculty of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi China ✉ Corresponding author. # Contributed equally. Received 2025 May 23; Accepted 2026 Feb 10; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ . PMC Copyright notice PMCID: PMC13065825  PMID: 41764165 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. Subject terms: Computational models, Data mining Spatial omics maps cell positions, but how they hierarchically assemble into functional tissues remains unclear. Here, the authors present HRCHY-CytoCommunity, a tool that decodes these multi-level organisations, linking them to biological mechanisms and clinical outcomes. 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-seq 4 , and Slide-seq 5 , 6 , provide transcriptome-wide gene expression profiling while preserving spatial context. Imaging-based spatial transcriptomics and proteomics technologies like MERFISH 7 , CODEX 8 , IMC 9 , and MIBI-TOF 10 , 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 a multi-level organization is therefore essential for understanding 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 a 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. Figure 1 illustrates the identification process of tissue structures at two hierarchical levels. Fig. 1. Overview of the HRCHY-CytoCommunity framework. Open in a new tab The framework formulates hierarchical tissue structure identification as a hierarchical community detection problem on a cell-cell proximity graph. HRCHY-CytoCommunity includes two modules: a a soft hierarchical structure assignment base module and b a robust hierarchical structure identification module. a The base module begins with a cell-cell proximity graph, where nodes represent cells, and their attributes are one-hot encoded cell-type vectors. A graph convolution layer followed by a fully-connected layer transforms node attribute vectors into soft fine-grained cellular neighborhood (CN) assignment vectors. A differentiable graph pooling layer is then used to generate a coarsene,d completed graph, where each pooled node represents a CN. Adaptive edge pruning is applied to mitigate over-smoothing. Another graph convolution and fully-connected layers process the pooled node embeddings to produce soft coarse-grained tissue compartment (TC) assignment vectors. The entire process is optimized by a graph minimum cut (MinCut)-based loss function, denoted as L Base . b To enhance the accuracy and robustness of hierarchical tissue structure identification, M perturbed cell-cell proximity graphs are generated via DropNode (random masking of node attribute vectors). Each graph is processed by the base module to produce soft fine-grained CN and coarse-grained TC assignment matrices. Consistency and balance regularization during training enforces stability across these M assignments. The total loss, denoted as L Total , combines the base MinCut loss ( L Base ) with a consistency regularization term ( L Consis , which minimizes divergence across M assignments) and a balance regularization term ( L Balance , which prevents cluster collapse). Finally, this module yields stable assignment matrices, from which final robust hierarchical tissue structures are derived. 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 the 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 a 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 are 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 36 . 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 37 – 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 ). Fig. 2. Performance evaluation of HRCHY-CytoCommunity on coarse-grained TC identification. Open in a new tab a Spatial maps of healthy mouse spleen samples generated by CODEX technology. Cells are colored by cell types (left) and manually annotated compartments (right), which serve as the ground truth for coarse-grained TCs. b Coarse-grained TCs identified by HRCHY-CytoCommunity (left) and NeST (right). c Coarse-grained TCs identified by non-hierarchical methods, including CellCharter, NicheCompass, GraphST, and SpaGCN. d Adjusted Mutual Information (AMI) and Macro-F1 scores calculated using manual compartment annotations. e Spatial maps of human CRC samples (8 μm bins) generated by Visium HD technology. Bins are colored by cell types (left) and manually annotated coarse-grained structures (right). f Coarse-grained TCs identified by HRCHY-CytoCommunity (left) and NeST (right). g Coarse-grained TCs identified by non-hierarchical methods, including CellCharter, NicheCompass, and SpaSEG. h AMI and Macro-F1 scores calculated using manual compartment annotations. In both d and h , each point corresponds to the performance on an individual sample, with horizontal bars indicating the mean performance across n = 3 samples. Points from the same sample are connected by grey dashed lines. P -values were calculated using one-sided paired t -tests. *, P -value < 0.05; **, P -value < 0.01; ***, P -value < 0.001. Unidentified indicates cells not assigned to any TC. Unmatched denotes TCs that could not be aligned with manual annotations. Source data are provided as a Source Data file. 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 splits 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 ). Fig. 3. Performance evaluation of HRCHY-CytoCommunity on fine-grained CN identification. Open in a new tab a Spatial maps of mouse hypothalamic preoptic region samples generated by MERFISH technology. Cells are colored by cell types (left) and manually annotated hypothalamic nuclei (right), which were generated based on the outlines of hypothalamic nuclei (middle) and serve as the ground truth for fine-grained CNs. b Fine-grained CNs identified by HRCHY-CytoCommunity (left) and NeST (right). c Fine-grained CNs identified by non-hierarchical methods, including CellCharter, NicheCompass, GraphST, SpaGCN, SpaSEG, and Giotto Suite. d Adjusted Mutual Information (AMI) and Macro-F1 scores calculated using manually annotated hypothalamic nuclei. Each point corresponds to the performance on an individual sample, with horizontal bars indicating the mean performance across n = 5 samples. Points from the same sample are connected by grey dashed lines. P -values were calculated using one-sided paired t -tests. *, P -value < 0.05; **, P -value < 0.01; ***, P -value < 0.001. Unidentified indicates cells not assigned to any CN. Source data are provided as a Source Data file. 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 its 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 the 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 biologically 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 , 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 (Fig. 4a, 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 (Fig. 4a, b , right panels). Fig. 4. Evaluation of HRCHY-CytoCommunity’ s cross-platform generalization using the mouse hippocampus Slide-seq V2 and human breast cancer Visium V1 datasets. Open in a new tab a Coarse-grained TC reference from the Allen Mouse Brain Atlas 42 (left), and TCs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the mouse hippocampus. b Fine-grained CN reference from the Allen Mouse Brain Atlas (left), and CNs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the mouse hippocampus. V3, third ventricle; MH, medial habenula; LH, lateral habenula; LP, lateral posterior nucleus of the thalamus; LD, lateral dorsal nucleus of the thalamus. c Manual annotation of coarse-grained TCs from Xu et al. 43 . (left), and TCs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the human breast cancer. d Manual annotation of fine-grained CNs from Xu et al. (left), and CNs identified by HRCHY-CytoCommunity (middle) and NeST (right) in the human breast cancer. e Heatmap showing cell-type enrichment scores for CNs identified by HRCHY-CytoCommunity. Enrichment score was defined as -log 10 (adjusted P -value). P -values were computed using one-sided Wilcoxon rank-sum tests and adjusted with the Benjamini-Hochberg method. *, P -value < 0.05; **, P -value < 0.01; ***, P -value < 0.001. Exact adjusted P -values are provided in the Source Data file. f Differential gene expression analysis of CN-8 (identified by HRCHY-CytoCommunity) versus other CNs. Each point represents a gene, with the y-axis showing -log 10 (adjusted P -value) and the x-axis showing log 2 (FoldChange). FC, Fold Change. P -values were computed using two-sided Wilcoxon rank-sum tests, with significance thresholds set at |log 2 FC | > 0.25 and P -value < 0.05. Unidentified indicates cells not assigned to any TC or CN. Unmatched denotes TCs or CNs that could not be aligned with manual annotations. Red arrowheads indicate the IDC_3 region manually annotated by Xu et al., which was reclassified by both HRCHY-CytoCommunity and NeST as a non-tumor CN. Source data are provided as a Source Data file. 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 (Fig. 4c, 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 IDC_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 JCHAIN ), 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 clear 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 ). Fig. 5. Cross-sample analysis of hierarchical tissue structures in the TNBC MIBI-TOF dataset. Open in a new tab a Representative single-cell spatial maps from four TNBC patients, showing compartmentalized tumors with clear spatial segregation between immune and neoplastic cells. Cells are colored by cell types (first column), coarse-grained TCs identified by HRCHY-CytoCommunity (second column) and NeST (third column), and fine-grained CNs identified by HRCHY-CytoCommunity (fourth column) and NeST (fifth column). Unidentified indicates cells not assigned to any TC or CN. b Proportion of neoplastic and immune cells correctly assigned to corresponding TCs. Each point corresponds to the performance on an individual sample, with horizontal bars indicating the mean performance across n = 15 samples. Points from the same sample are connected by grey dashed lines. P -value was calculated using a one-sided paired t -test. c Unified CN maps for patient 4 and patient 9. Cells are colored by unified CNs. Purple dashed lines outline neoplastic cell-dominated TCs, while remaining regions denote immune cell-dominated TCs. d Heatmaps showing cell-type enrichment scores for each unified CN in patient 4 and patient 9. Enrichment score was defined as -log 10 (adjusted P -value). P -values were computed using hypergeometric tests for enrichment (over-representation) and adjusted with the Benjamini-Hochberg method. Exact adjusted P -values are provided in the Source Data file. e Hierarchical clustering trees based on unified CN composition profiles (left) and cell-type composition profiles (right). Source data are provided as a Source Data file. 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 (Fig. 5c, 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 ). TC-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 9 . Fig. 6. Hierarchical tissue structure-based stratification of breast cancer patients. Open in a new tab a Heatmaps showing the fraction of patients in each group containing each unified TC (top), and average cell-type enrichment scores for each unified TC (bottom). b Kaplan-Meier survival curves of 263 breast cancer patients, categorized into five groups based on unified TCs. c Forest plot of hazard ratios (HRs) and 95% confidence intervals (CIs) from multivariable Cox proportional hazards regression (adjusting for age, tumor grade, and tumor size) evaluating the association between unified TC-13 and overall survival ( N = 263 patients, 73 events). Presence of unified TC-13 was significantly associated with improved survival (HR = 0.20, P -value = 0.031). d Forest plot of HRs and 95% CIs from multivariable Cox regression for unified CN-21 within Group 1 ( N = 122 patients, 34 events). Presence of unified CN-21 was significantly associated with worse survival (HR = 2.2, P -value = 0.026). e Forest plot of HRs and 95% CIs from multivariable Cox regression for unified CN-34 within Group 1 ( N = 122 patients, 34 events). Presence of unified CN-34 was significantly associated with worse survival (HR = 5.4, P -value < 0.001). In ( c – e ), square markers indicate HR point estimates, and horizontal error bars represent the corresponding 95% CIs. f Heatmaps showing cell-type enrichment scores for unified CN-21 and CN-34 across patients. Significantly enriched cell types in these two CNs are highlighted by red dashed boxes. In both a and f enrichment score was defined as -log 10 (adjusted P -value). P -values were computed using hypergeometric tests for enrichment (over-representation) and adjusted with the Benjamini-Hochberg method. *, P -value < 0.05; **, P -value < 0.01; ***, P -value < 0.001. Exact adjusted P -values are provided in the Source Data file. g Representative single-cell spatial maps from patients containing unified CN-21. Cells are colored by cell types (top) and unified CNs (bottom). h Representative single-cell spatial maps from patients containing unified CN-34. Cells are colored by cell types (top) and unified CNs (bottom). Source data are provided as a Source Data file. 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 (Fig. 6d, 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 Fibronectin hi , Vimentin hi , and SMA hi Vimentin hi 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. Discussion We present HRCHY-CytoCommunity, a GNN-based framework designed to identify hierarchical tissue structures based on cell types and their spatial locations. By leveraging differentiable graph pooling, adaptive edge pruning, and consistency and balance regularization during training, HRCHY-CytoCommunity provides a robust, end-to-end approach that directly links single-cell features to multi-scale tissue architecture. Unlike existing hierarchical methods such as NeST 34 , which relies on gene expression hotspots, HRCHY-CytoCommunity uses cell-type annotations as core features. This design makes it widely applicable to diverse spatial omics technologies, especially those with limited gene or protein features, that can characterize cellular phenotypes from complementary perspectives. Furthermore, in contrast to NeST, HRCHY-CytoCommunity produces fully nested hierarchical assignments with complete cellular coverage, offering clearer biological interpretability and enabling a more comprehensive understanding of how individual cells collectively build tissues and organs. Different from non-hierarchical methods that rely on post hoc adjustment of clustering resolution to infer pseudo-hierarchies, HRCHY-CytoCommunity integrates the learning of multiple organizational levels within a unified end-to-end framework. This joint optimization approach enables information to flow bidirectionally across scales during training, promoting the identification of hierarchically consistent and functionally relevant tissue architectures. Extensive benchmarking on spatial omics datasets from various tissues, technologies, and modalities demonstrates that HRCHY-CytoCommunity consistently outperforms NeST and other non-hierarchical methods in terms of both accuracy and biological relevance. The incorporation of consistency and balance regularization further enhances its robustness and reproducibility. Moreover, through efficient sparse graph operations, HRCHY-CytoCommunity achieves high computational scalability, making it well-suited for large-scale and high-resolution spatial datasets. To support cross-sample integration analysis, HRCHY-CytoCommunity includes an additional module that performs cell-type enrichment-based clustering, generating a unified set of nested multicellular structures across all samples. This addresses the challenge of inconsistent structural labels between samples. Applied to a TNBC dataset, the method successfully identified both conserved and patient-specific tissue architectures. Furthermore, using hierarchical tissue structures as features in survival analysis, we hierarchically stratified breast cancer patients from an IMC dataset into prognostically distinct groups, underscoring the clinical relevance of multi-scale spatial organization and suggesting new avenues for exploring hierarchical tumor microenvironments and related immunotherapeutic strategies. As HRCHY-CytoCommunity explicitly models tissue organization as a hierarchy, it is particularly well-suited for organs with clearly stratified anatomical or functional layers. Its performance may be limited in tissues lacking strong hierarchical organization, such as liver, where metabolic zonation is continuously distributed along the porto-central axis 52 . Moreover, while the use of discrete cell-type annotations enhances the method’s generalizability across platforms, it may not fully capture structures defined by continuous molecular gradients, transient states, or cell-cell communication signals. Future versions could incorporate multi-view learning to integrate such continuous features. In the meantime, methods like NicheCompass 20 (focused on cellular communication-driven structures) and ONTraC 53 (designed for spatially continuous niche trajectories) can complement HRCHY-CytoCommunity to provide a more comprehensive perspective on tissue organization. Further extensions could also incorporate graph sampling or alignment techniques to improve cross-sample integration within an end-to-end analytical workflow. Lastly, although this study focused on a two-level hierarchy, the underlying framework of HRCHY-CytoCommunity can be readily extended to deeper hierarchies (Supplementary Note 8 ). Future work may explore adaptive mechanisms to automatically infer the optimal hierarchical depth or to model tissues with non-uniform organizational principles. As spatial omics technologies continue to evolve, there is a growing need for versatile and scalable computational methods that can decipher tissue organization across platforms with varying resolution and feature throughput. By focusing on cell-type spatial maps, a common output across spatial omics technologies, HRCHY-CytoCommunity represents a general and scalable framework for simultaneously identifying multi-level hierarchical tissue structures. It provides a foundational tool for decoding principles of tissue organization by systematically bridging spatial hierarchies from single cells to multicellular modules, and ultimately to intact tissues and organs. 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 Z ( 1 ) ∈ R n × d , formulated as Z ( 1 ) = ReLU GNN 1 X ( 1 ) , A ( 1 ) ; θ GNN 1 1 where each row of Z ( 1 ) is a learned d -dimensional embedding vector of a cell. X ( 1 ) ∈ R n × t denotes the initial cell-type attribute matrix, and t is the total number of cell types. A ( 1 ) ∈ { 0, 1 } n × n is the adjacency matrix of the undirected KNN graph. The graph convolution operator was defined as z i ( 1 ) = ReLU θ 1 x i ( 1 ) + θ 2 ∑ j ∈ N ( i ) x j ( 1 ) 2 where x i is the embedding vector of node i , and N ( i ) denotes the first-order neighborhood derived from the matrix A ( 1 ) . θ 1 and θ 2 are trainable parameters in the graph neural network GNN 1 . The embedding dimension d was empirically set to 128 in this study. 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 Z ( 1 ) ∈ R n × d into a soft fine-grained tissue structure assignment matrix S ( 1 ) ∈ R n × c 1 , which was formulated as below. S ( 1 ) = Softmax FC 1 Z ( 1 ) ; θ FC 1 3 where each element of S ( 1 ) represents the probability of a cell (row) belonging to one of the c 1 fine-grained tissue structures (column). θ FC 1 represents trainable parameters in the fully-connected neural network FC 1 . The hyperparameter c 1 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 c 1 ). To further identify coarse-grained tissue structures, we applied a differentiable graph pooling layer 35 , 55 to generate a coarsened graph with c 1 pooled nodes, where each pooled node represents a fine-grained tissue structure. The adjacency matrix A ( 2 ) ∈ R c 1 × c 1 and node embedding matrix X ( 2 ) ∈ R c 1 × d of the coarsened graph were formulated as follows. A ( 2 ) = S ( 1 ) T A ( 1 ) S ( 1 ) 4 X ( 2 ) = S ( 1 ) T Z ( 1 ) 5 The resulting coarsened graph is fully-connected with self-loops. Self-loops were removed and the adjacency matrix was normalized as A ^ = A ( 2 ) − I c 1 d i a g ( A ( 2 ) ) 6 A p o o l = D ^ − 1 2 A ^ D ^ − 1 2 7 where D ^ is the degree matrix of A ^ . Each entry in A p o o l 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 T e d g e _ p r u n i n g . Only edges with weights in the adjacency matrix A p o o l greater than this threshold are retained and their weights are reset to 1. A ^ i j ( 2 ) = 1 i f A i j p o o l > T e d g e _ p r u n i n g 0 o t h e r w i s e 8 Since excessive or insufficient edge pruning may affect model performance (Supplementary Note 3 ; Supplementary Figs. 13b – 15b ), the adaptive threshold was defined as T e d g e _ p r u n i n g = 1 c 1 − 1 9 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: Z ( 2 ) = ReLU GNN 2 X ( 2 ) , A ^ ( 2 ) ; θ GNN 2 10 S ( 2 ) = Softmax FC 2 Z ( 2 ) ; θ FC 2 11 Each element of S ( 2 ) ∈ R c 1 × c 2 represents the probability of a fine-grained structure (row) belonging to one of the c 2 coarse-grained structures (column). Z ( 2 ) ∈ R c 1 × d is the updated embedding matrix of pooled nodes in the coarsened graph. θ GNN 2 and θ FC 2 represent trainable parameters in the graph neural network GNN 2 and fully-connected neural network FC 2 , respectively. The hyperparameter c 2 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 c 2 ). The loss function with regard to the soft hierarchical tissue structure assignment (base module) was defined as follows. L B a s e = α × L F i n e + 1 − α × L C o a r s e 12 where α is a weight parameter used for balancing the fine-grained loss L F i n e and the coarse-grained loss L C o a r s e . 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 L F i n e and L C o a r s e use the graph MinCut-based formulation, optimizing the matrix S ( 1 ) and S ( 2 ) , respectively. L F i n e = − ∑ j = 1 c 1 S 1 T A 1 S 1 j j ∑ j = 1 c 1 S 1 T D 1 S 1 j j + S ( 1 ) T S ( 1 ) ∣ ∣ S ( 1 ) T S ( 1 ) ∣ ∣ F − I c 1 c 1 F 13 L C o a r s e = − ∑ k = 1 c 2 S 2 T A 2 S 2 k k ∑ k = 1 c 2 S 2 T D 2 S 2 k k + S ( 2 ) T S ( 2 ) ∣ ∣ S ( 2 ) T S ( 2 ) ∣ ∣ F − I c 2 c 2 F 14 where D ( 1 ) ∈ R n × n and D ( 2 ) ∈ R c 1 × c 1 are degree matrices derived from the undirected KNN graph and the coarsened graph, respectively. Both L F i n e and L C o a r s e consist of two terms. The first term encourages clustering of strongly connected nodes, and the second term encourages orthogonal and balanced tissue structure assignments 35 . 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 A ( 1 ) and node attribute matrix X ( 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 X ( 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 i , a binary mask ϵ i was sampled from a Bernoulli distribution ϵ i ~ Bernoulli ( 1 − δ ) , where δ is a predefined DropNode rate. The perturbed feature vector was then computed as x ~ i = ϵ i ⋅ x i , where x i represents the original feature vector of node i . The perturbed feature matrix X ~ was used only during training. During inference, features were rescaled by multiplying ( 1 − δ ) to match the expected activation magnitude of X ~ . 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: s ¯ i = 1 m ∑ m = 1 M s i ( m ) 15 where s i ( m ) denotes the predicted tissue structure assignment distribution of node i in the m -th perturbed graph. Finally, the consistency regularization loss was defined as the average L 2 distance between the assignment distributions from each perturbed graph and the distribution center: L C o n s i s = 1 M ∑ m = 1 M ∑ i = 1 n s i ( m ) − s ¯ i 2 2 16 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 a collapse, we introduced an entropy-based balance regularization term: L B a l a n c e = 1 − 1 M ∑ m = 1 M H m ( 2 ) log c 2 17 H m ( 2 ) = − ∑ k = 1 c 2 s m , k ( 2 ) log s m , k ( 2 ) 18 where H m ( 2 ) is the entropy of the cluster assignment distribution for the m -th perturbation, and c 2 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 L T o t a l = L B a s e + λ 1 × L C o n s i s + λ 2 × L B a l a n c e 19 where L B a s e is the base loss from the soft hierarchical tissue structure assignment, L C o n s i s is the consistency regularization loss, and L B a l a n c e is the entropy-based balance regularization loss. The hyperparameters λ 1 and λ 2 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: F M I = T P ( T P + F P ) ( T P + F N ) 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 the 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 12 ( 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 a 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 36 . For the mouse hypothalamic preoptic region MERFISH dataset, the outlines of hypothalamic nuclei regions were obtained from the original study 7 , and the ground-truth assignments of cells to CNs were manually annotated based on these outlines. For the mouse ICH Stereo-seq dataset, the CN annotations were obtained from the original study 60 . 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. F 1 s c o r e = 2 × ( P r e c i s i o n × R e c a l l ) P r e c i s i o n + R e c a l l 21 P r e c i s i o n = T P T P + F P 22 R e c a l l = T P T P + F N 23 A M I = I G T ; T S − E { I G T ; T S } 1 2 H G T + H T S − E { I G T ; T S } 24 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. E { I G T ; T S } denotes the expected mutual information between the ground-truth assignments of cells (GT) and the tissue structure predictions (TS). H G T and H ( T S ) 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 61 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 -log 10 ( 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 62 . Survival analysis To validate the 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 9 . 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. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Supplementary information Supplementary Information (17.4MB, pdf) 41467_2026_70069_MOESM2_ESM.pdf (27.5KB, pdf) Description of Additional Supplementary Files Supplementary Data 1 (15.2KB, xlsx) Reporting Summary (86.4KB, pdf) Transparent Peer Review file (54.6MB, pdf) Source data Source Data (2MB, xlsx) Acknowledgements This work was supported by a National Natural Science Foundation of China (NSFC) grant No. 62422211 and a Scientific Research Innovation Capability Support Project for Young Faculty No. SRICSPYF-ZY2025003 to Y. H., NSFC grants no. 62550005, 62132015, and No. U22A2037 to L. G., and an NSFC grant no. 62302386 and a National Key Research and Development Program of China No. 2024YFC2707102 to J. Lin. We thank the Key Laboratory of Computational Bioinformatics of Xi’an at Xidian University for providing computing support. Author contributions Y. H. and L. G. conceived and designed the study. R. X., Z. W., and Y. H. designed the HRCHY-CytoCommunity algorithm. R. X. and Z. W. implemented the HRCHY-CytoCommunity algorithm. J. Liu, H. X., and Y. X. provided additional input during the method development. R. X., Z. W., J. Liu, and Y. H. performed the data analysis. R. X. and Z. W. provided support for the software package development. Y. H., L. G., and J. Lin supervised the study. R. X., Y. H., Z. W., J. Lin, and L. G. wrote the manuscript. Peer review Peer review information Nature Communications thanks Sergio Marco Salas, who co-reviewed with Soroor Hediyeh-zadeh, Jun Ding, and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available. Data availability All data used in this study are publicly available. This study used the following eight publicly available spatial omics datasets, including a mouse spleen CODEX dataset ( https://data.mendeley.com/datasets/zjnpwh8m5b/1 ), a mouse hypothalamic preoptic region MERFISH dataset ( https://datadryad.org/stash/dataset/doi:10.5061/dryad.8t8s248 ), a human triple-negative breast cancer MIBI-TOF dataset ( https://mibi-share.ionpath.com ), a human breast cancer IMC dataset ( https://zenodo.org/record/3518284#.Y2UQ0-xBybg ), a human CRC Visium HD dataset ( https://www.10xgenomics.com/platforms/visium/product-family/dataset-human-crc ), a human breast cancer Visium V1 dataset ( https://www.10xgenomics.com/datasets/human-breast-cancer-block-a-section-1-1-standard-1-1-0 ), a mouse hippocampus Slide-seq V2 dataset ( https://singlecell.broadinstitute.org/single_cell/study/SCP815/highly-sensitive-spatial-transcriptomics-at-near-cellular-resolution-with-slide-seqv2 ), and a mouse ICH Stereo-seq dataset ( https://db.cngb.org/stomics/stmich/ ). Two single-cell transcriptomic datasets were used as references for cell-type deconvolution in the spot-based spatial transcriptomics data, including a human breast cancer scRNA-seq dataset (GEO: GSE176078 ; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE176078 ), and a mouse hippocampus scRNA-seq dataset ( https://singlecell.broadinstitute.org/single_cell/study/SCP948/robust-decomposition-of-cell-type-mixtures-in-spatial-transcriptomics#study-download ). No new data were generated in this study. Source data are provided with this paper. Code availability HRCHY-CytoCommunity is an open-access Python package available in the GitHub repository at https://github.com/huBioinfo/HRCHY-CytoCommunity , under the MIT license. The specific version of the code associated with this publication is archived in Zenodo and is accessible via https://zenodo.org/records/18137898 63 . Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. These authors contributed equally: Runzhi Xie, Zekun Wang. Contributor Information Jiadong Lin, Email: [email protected]. Yuxuan Hu, Email: [email protected]. Lin Gao, Email: [email protected]. Supplementary information The online version contains supplementary material available at 10.1038/s41467-026-70069-z. References 1. Bressan, D., Battistoni, G. & Hannon, G. J. The dawn of spatial omics. Science 381 , eabq4964 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 2. Ståhl, P. L. et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science 353 , 6294 (2016). [ DOI ] [ PubMed ] [ Google Scholar ] 3. Oliveira, M. F. D. et al. High-definition spatial transcriptomic profiling of immune cell populations in colorectal cancer. Nat. Genet. 57 , 1512–1523 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 4. Chen, A. et al. Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays. Cell 185 , 1777–1792.e21 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 5. Rodriques, S. G. et al. Slide-seq: A scalable technology for measuring genome-wide expression at high spatial resolution. Science 363 , 6434 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 6. Stickels, R. R. et al. Highly sensitive spatial transcriptomics at near-cellular resolution with Slide-seqV2. Nat. Biotechnol. 39 , 313–319 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 7. Moffitt, J. R. et al. Molecular, spatial, and functional single-cell profiling of the hypothalamic preoptic region. Science 362 , eaau5324 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 8. Schürch, C. M. et al. Coordinated cellular neighborhoods orchestrate antitumoral immunity at the colorectal cancer invasive front. Cell 182 , 1341–1359.e19 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 9. Jackson, H. W. et al. The single-cell pathology landscape of breast cancer. Nature 578 , 615–620 (2020). [ DOI ] [ PubMed ] [ Google Scholar ] 10. Keren, L. et al. A structured tumor-immune microenvironment in triple negative breast cancer revealed by multiplexed ion beam imaging. Cell 174 , 1373–1387.e19 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Rao, A., Barkley, D., França, G. S. & Yanai, I. Exploring tissue architecture using spatial transcriptomics. Nature 596 , 211–220 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 12. Zhu, Q., Shah, S., Dries, R., Cai, L. & Yuan, G.-C. Identification of spatially associated subpopulations by combining scRNAseq and sequential fluorescence in situ hybridization data. Nat. Biotechnol. 36 , 1183–1190 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 13. Li, J., Chen, S., Pan, X., Yuan, Y. & Shen, H.-B. Cell clustering for spatial transcriptomics data with graph neural networks. Nat. Comput Sci. 2 , 399–408 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 14. Dong, K. & Zhang, S. Deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder. Nat. Commun. 13 , 1739 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 15. Long, Y. et al. Spatially informed clustering, integration, and deconvolution of spatial transcriptomics with GraphST. Nat. Commun. 14 , 1155 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 16. Hu, Y. et al. Unsupervised and supervised discovery of tissue cellular neighborhoods from cell phenotypes. Nat. Methods 21 , 267–278 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 17. Ma, Y. & Zhou, X. Accurate and efficient integrative reference-informed spatial domain detection for spatial transcriptomics. Nat. Methods 21 , 1231–1244 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 18. Long, Y. et al. Deciphering spatial domains from spatial multi-omics with SpatialGlue. Nat. Methods 21 , 1658–1667 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 19. Singhal, V. et al. BANKSY unifies cell typing and tissue domain segmentation for scalable spatial omics data analysis. Nat. Genet. 56 , 431–441 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 20. Birk, S. et al. Quantitative characterization of cell niches in spatially resolved omics data. Nat. Genet. 57 , 897–909 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 21. Varrone, M., Tavernari, D., Santamaria-Martínez, A., Walsh, L. A. & Ciriello, G. CellCharter reveals spatial cell niches associated with tissue remodeling and cell plasticity. Nat. Genet. 56 , 74–84 (2024). [ DOI ] [ PubMed ] [ Google Scholar ] 22. Hu, J. et al. SpaGCN: Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network. Nat. Methods 18 , 1342–1351 (2021). [ DOI ] [ PubMed ] [ Google Scholar ] 23. Bai, Y. et al. SpaSEG: unsupervised deep learning for multi-task analysis of spatially resolved transcriptomics. Genome Biol. 26 , 230 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 24. Chen, J. G. et al. Giotto Suite: a multiscale and technology-agnostic spatial multiomics analysis ecosystem. Nat. Methods 22 , 2052–2064 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Bhate, S. S., Barlow, G. L., Schürch, C. M. & Nolan, G. P. Tissue schematics map the specialization of immune tissue motifs and their appropriation by tumors. Cell Syst. 13 , 109–130.e6 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Hickey, J. W. et al. Organization of the human intestine at single-cell resolution. Nature 619 , 572–584 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 27. Nettekoven, C. et al. A hierarchical atlas of the human cerebellum for functional precision mapping. Nat. Commun. 15 , 8376 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 28. de Vries, H.-M., Ottenhof, S. R., Horenblas, S., van der Heijden, M. S. & Jordanova, E. S. Defining the tumor microenvironment of penile cancer by means of the cancer immunogram. Eur. Urol. Focus 5 , 718–721 (2019). [ DOI ] [ PubMed ] [ Google Scholar ] 29. Li, J. et al. Remodeling of the immune and stromal cell compartment by PD-1 blockade in mismatch repair-deficient colorectal cancer. Cancer Cell 41 , 1152–1169.e7 (2023). [ DOI ] [ PubMed ] [ Google Scholar ] 30. Fridman, W. H. et al. B cells and tertiary lymphoid structures as determinants of tumour immune contexture and clinical outcome. Nat. Rev. Clin. Oncol. 19 , 441–457 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 31. Schumacher, T. N. & Thommen, D. S. Tertiary lymphoid structures in cancer. Science 375 , eabf9419 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 32. Teillaud, J.-L., Houel, A., Panouillot, M., Riffard, C. & Dieu-Nosjean, M.-C. Tertiary lymphoid structures in anticancer immunity. Nat. Rev. Cancer 24 , 629–646 (2024). [ DOI ] [ PubMed ] [ Google Scholar ] 33. Sato, Y., Silina, K., van den Broek, M., Hirahara, K. & Yanagita, M. The roles of tertiary lymphoid structures in chronic diseases. Nat. Rev. Nephrol. 19 , 525–537 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 34. Walker, B. L. & Nie, Q. NeST: nested hierarchical structure identification in spatial transcriptomic data. Nat. Commun. 14 , 6554 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 35. Bianchi, F. M., Grattarola, D. & Alippi, C. Spectral Clustering with Graph Neural Networks for Graph Pooling. in Proceedings of the 37th International Conference on Machine Learning 874–883 (PMLR, 2020). 36. Goltsev, Y. et al. Deep profiling of mouse splenic architecture with CODEX Multiplexed Imaging. Cell 174 , 968–981.e15 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 37. Gonzalez, M., Mackay, F., Browning, J. L., Kosco-Vilbois, M. H. & Noelle, R. J. The Sequential Role of Lymphotoxin and B Cells in the Development of Splenic Follicles. J. Exp. Med. 187 , 997–1007 (1998). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 38. Mebius, R. E. & Kraal, G. Structure and function of the spleen. Nat. Rev. Immunol. 5 , 606–616 (2005). [ DOI ] [ PubMed ] [ Google Scholar ] 39. Cerutti, A., Cols, M. & Puga, I. Marginal zone B cells: virtues of innate-like antibody-producing lymphocytes. Nat. Rev. Immunol. 13 , 118–132 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 40. Paxinos and Franklin’s the Mouse Brain in Stereotaxic Coordinates . (Academic Press, 2019). 41. Cable, D. M. et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat. Biotechnol. 40 , 517–526 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 42. Wang, Q. et al. The Allen Mouse Brain Common Coordinate Framework: A 3D Reference Atlas. Cell 181 , 936–953.e20 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 43. Xu, H. et al. Unsupervised spatially embedded deep representation of spatial transcriptomics. Genome Med. 16 , 12 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 44. Hollern, D. P. et al. B Cells and T follicular helper cells mediate response to checkpoint inhibitors in high mutation burden mouse models of breast cancer. Cell 179 , 1191–1206.e21 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 45. Gentles, A. J. et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat. Med 21 , 938–945 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 46. Meylan, M. et al. Tertiary lymphoid structures generate and propagate anti-tumor antibody-producing plasma cells in renal cell cancer. Immunity 55 , 527–541.e5 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 47. Chhabra, Y. & Weeraratna, A. T. Fibroblasts in cancer: Unity in heterogeneity. Cell 186 , 1580–1609 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 48. Liu, Y. et al. Conserved spatial subtypes and cellular neighborhoods of cancer-associated fibroblasts revealed by single-cell spatial multi-omics. Cancer Cell 43 , 905–924.e6 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 49. Ye, F. et al. Cancer-associated fibroblasts facilitate breast cancer progression through exosomal circTBPL1-mediated intercellular communication. Cell Death Dis. 14 , 471 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 50. Barbazan, J. et al. Cancer-associated fibroblasts actively compress cancer cells and modulate mechanotransduction. Nat. Commun. 14 , 6966 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 51. Roerden, M. & Spranger, S. Cancer immune evasion, immunoediting and intratumour heterogeneity. Nat. Rev. Immunol. 25 , 353–369 (2025). [ DOI ] [ PubMed ] [ Google Scholar ] 52. Kietzmann, T. Metabolic zonation of the liver: The oxygen gradient revisited. Redox Biol. 11 , 622–630 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 53. Wang, W., Zheng, S., Shin, S. C., Chávez-Fuentes, J. C. & Yuan, G.-C. ONTraC characterizes spatially continuous variations of tissue microenvironment through niche trajectory analysis. Genome Biol. 26 , 117 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 54. Morris, C. et al. Weisfeiler and Leman go neural: higher-order graph neural networks. Proc. AAAI Conf. Artif. Intell. 33 , 4602–4609 (2019). [ Google Scholar ] 55. Ying, Z. et al. Hierarchical graph representation learning with differentiable pooling. Adv. Neural Inf. Process. Syst. 31 (2018). 56. Feng, W. et al. Graph random neural networks for semi-supervised learning on graphs. Adv. neural Inf. Process. Syst. 33 , 22092–22103 (2020). [ Google Scholar ] 57. Fowlkes, E. B. & Mallows, C. L. A Method for Comparing Two Hierarchical Clusterings. J. Am. Stat. Assoc. 78 , 553–569 (1983). [ Google Scholar ] 58. Lotfollahi, M. et al. Mapping single-cell data to reference atlases by transfer learning. Nat. Biotechnol. 40 , 121–130 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 59. Lopez, R., Regier, J., Cole, M. B., Jordan, M. I. & Yosef, N. Deep generative modeling for single-cell transcriptomics. Nat. Methods 15 , 1053–1058 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 60. Xiang, R. et al. Spatiotemporal transcriptomic maps of mouse intracerebral hemorrhage at single-cell resolution. Neuron (2025). [ DOI ] [ PubMed ] 61. Romano, S., Vinh, N. X., Bailey, J. & Verspoor, K. Adjusting for Chance Clustering Comparison Measures. J. Mach. Learn. Res. 17 , 1–32 (2016). [ Google Scholar ] 62. Benjamini, Y. & Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc.: Ser. B (Methodol.) 57 , 289–300 (1995). [ Google Scholar ] 63. Xie, R. et al. HRCHY-CytoCommunity. Zenodo 10.5281/zenodo.18137898 (2026). Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplementary Information (17.4MB, pdf) 41467_2026_70069_MOESM2_ESM.pdf (27.5KB, pdf) Description of Additional Supplementary Files Supplementary Data 1 (15.2KB, xlsx) Reporting Summary (86.4KB, pdf) Transparent Peer Review file (54.6MB, pdf) Source Data (2MB, xlsx) Data Availability Statement All data used in this study are publicly available. This study used the following eight publicly available spatial omics datasets, including a mouse spleen CODEX dataset ( https://data.mendeley.com/datasets/zjnpwh8m5b/1 ), a mouse hypothalamic preoptic region MERFISH dataset ( https://datadryad.org/stash/dataset/doi:10.5061/dryad.8t8s248 ), a human triple-negative breast cancer MIBI-TOF dataset ( https://mibi-share.ionpath.com ), a human breast cancer IMC dataset ( https://zenodo.org/record/3518284#.Y2UQ0-xBybg ), a human CRC Visium HD dataset ( https://www.10xgenomics.com/platforms/visium/product-family/dataset-human-crc ), a human breast cancer Visium V1 dataset ( https://www.10xgenomics.com/datasets/human-breast-cancer-block-a-section-1-1-standard-1-1-0 ), a mouse hippocampus Slide-seq V2 dataset ( https://singlecell.broadinstitute.org/single_cell/study/SCP815/highly-sensitive-spatial-transcriptomics-at-near-cellular-resolution-with-slide-seqv2 ), and a mouse ICH Stereo-seq dataset ( https://db.cngb.org/stomics/stmich/ ). Two single-cell transcriptomic datasets were used as references for cell-type deconvolution in the spot-based spatial transcriptomics data, including a human breast cancer scRNA-seq dataset (GEO: GSE176078 ; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE176078 ), and a mouse hippocampus scRNA-seq dataset ( https://singlecell.broadinstitute.org/single_cell/study/SCP948/robust-decomposition-of-cell-type-mixtures-in-spatial-transcriptomics#study-download ). No new data were generated in this study. Source data are provided with this paper. HRCHY-CytoCommunity is an open-access Python package available in the GitHub repository at https://github.com/huBioinfo/HRCHY-CytoCommunity , under the MIT license. The specific version of the code associated with this publication is archived in Zenodo and is accessible via https://zenodo.org/records/18137898 63 . Articles from Nature Communications are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (5.3 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top