Tumor cell villages define the co-dependency of tumor and microenvironment in liver cancer

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

Abstract
Spatial cellular context is crucial in shaping intratumor heterogeneity. However, understanding how each tumor establishes its unique spatial landscape and what factors drive the landscape for tumor fitness remain significantly challenging. Here, we analyze over 2 million cells from 50 tumor biospecimens using spatial single-cell imaging and single-cell RNA sequencing. We develop a deep learning-based strategy to spatially map tumor cell states and their surrounding environmental architecture, and find that different tumor cell states can be organized into distinct clusters, or “villages,” each supported by unique microenvironments. Notably, tumor cell villages exhibit village-specific molecular co-dependencies between tumor cells and their microenvironment and are associated with patient outcomes. Perturbation of molecular co-dependencies via random spatial shuffling of the microenvironment results in destabilization of the corresponding villages. We validate our findings using single-cell, spatial, and bulk transcriptome data from 740 liver cancer patients. This study provides insights into understanding tumor spatial landscape and its impact on tumor aggressiveness.

Introduction
Intratumor heterogeneity (ITH) is a major hurdle for effective interventions across a broad range of cancers. It has been increasingly studied by single-cell transcriptome analysis, where each tumor may contain various transcriptomic states distinguished by unique gene expression patterns1–3. In pan-cancer single-cell transcriptome studies, diverse tumor cell transcriptomic states have been defined, including cell cycle, epithelial-to-mesenchymal transition (EMT), and major histocompatibility complex (MHC) II-related states, among others4,5. While these features of ITH are well observed, the mechanisms by which diverse tumor cells spatially coordinate and how each tumor establishes its unique spatial landscape remain poorly understood.

Spatial transcriptomic profiling of various cancers at the spot resolution (a mixture of cells in each spot) has begun to unravel the intratumor spatial landscape, where a non-random spatial distribution of diverse transcriptome-based tumor clusters is observed within individual tumors6–13. In ovarian cancer, a specific tumor cell state was found to spatially colocalize with tumor-infiltrating lymphocytes based on single-cell spatial transcriptome analysis14. These observations suggest that the spatial cellular context is crucial in shaping ITH and further motivated us to understand the spatial landscape of individual tumors.

A critical challenge is how to spatially define tumor cell landscapes and identify the drivers underlying diverse spatial architectures. Tumors from different patients evolve independently, leading to unique gene expression profiles that reflect patient-specific characteristics. Directly clustering spatially resolved tumor spots or cells within each tumor tends to capture patient-specific patterns, making it difficult to compare the landscapes across patients. Here, we introduce the concept of “villages” to define the distinct landscapes formed by tumor cells, where each village arises from the spatial coordination of diverse tumor cell transcriptomic states supported by specific local environments, termed the Spatial Dynamics Network (SDN), surrounding individual tumor cells. Tumor cell villages represent the smallest community within a tumor. This concept draws inspiration from human villages, which are small, close-knit communities where individuals share resources and responsibilities. Like human villages, tumor cell “villages” may use communal mechanisms to promote growth and enhance defense, reducing individual cell vulnerability. The introduction of this village concept may provide insights into how tumor cells establish distinct spatial landscapes within individual patients, while also enabling comparisons of these landscapes across patients. More importantly, the defined tumor cell villages may uncover unique molecular co-dependencies among cells that are crucial for orchestrating their organization within each village, offering opportunities for therapeutic interventions.

In this study, we focus on primary liver cancer (PLC) due to its well-known, extensive ITH, making it an ideal model to investigate the spatial relationships of diverse tumor cell populations within individual tumors1,6,15–18. We perform transcriptomic profiling of liver tumors at single-cell spatial resolution using CosMx™ SMI 1000-plex in situ spatial imaging platform19. Additionally, we conduct single-cell RNA sequencing (scRNA-seq) on tumors from the same patients to generate more comprehensive expression profiles of individual cells. In total, we profile ~2.4 million cells from 50 tumor biospecimens of 7 patients with PLC, including 4 with hepatocellular carcinoma (HCC) and 3 with intrahepatic cholangiocarcinoma (iCCA). We define tumor cell transcriptomic states and the SDNs surrounding individual tumor cells. We observe distinct spatial preferences of tumor cell states and the coordination of diverse tumor cells as tumor cell “villages” supported by unique SDNs. Tumor cell villages are found to be linked to patient outcomes. Notably, village-specific molecular co-dependencies between tumor cells and their microenvironment are observed, the perturbation of which leads to the destabilization of the corresponding tumor cell villages. We validate our findings using single-cell transcriptome data, 10× Visium spatial transcriptome data, and bulk transcriptome data from 766 samples of 740 PLC patients. This study may provide crucial insights into understanding diverse tumor spatial architectures and their impact on tumor aggressiveness.

Results
Single-cell spatial characterization of tumors from primary liver cancer patients
To comprehensively capture ITH and investigate its stability within individual tumors, we designed the study to collect samples from multiple locations within each tumor from PLC patients (Fig. 1A). In total, 50 samples were collected from the tumor core, tumor border, and adjacent non-tumor tissue of 7 PLC patients (4 HCC, 3 iCCA; Supplementary Data 1). Among them, 15 samples (whole tissue section) were applied for single-cell spatial transcriptome profiling using the CosMx™ SMI 1000-plex in situ spatial imaging platform19. The remaining samples were profiled by single-cell RNA sequencing (scRNA-seq) on tumors from the same patients to generate more comprehensive expression profiles of individual cells in our previous study using the 10× Genomics droplet-based 5′ assay20. In this study, we refer to the CosMx™ SMI data as “single-cell spatial data” and the scRNA-seq data as “single-cell data”. After quality control (Methods), spatial transcriptomic profiles of 2,347,589 cells were derived (Fig. 1B, Supplementary Fig. 1A, B). In addition, scRNA-seq data of 112,506 cells were generated (Fig. 1B)20. When mapping individual cells to their spatial locations, we observed a strong concordance among cell type annotations, the expression of cell type-specific marker genes, and protein staining, demonstrating a good representation of the tumor landscape with the single-cell spatial data (Fig. 1C–E, Supplementary Fig. 1B, C; Supplementary Data 2).

To further determine the quality of the single-cell spatial transcriptome data, we compared it with the scRNA-seq data from the same patients. We found strong correlations in gene expression between the two platforms, confirming a high quality of the spatial data (Supplementary Fig. 1D). Interestingly, we observed a striking difference in cell type composition from the two platforms. In the scRNA-seq data, T cells accounted for over 80% of the entire cell population (Fig. 1F). By contrast, epithelial cells represented the major cell type in the single-cell spatial data, which was consistent with the histological images of the tumors20. This discrepancy may be due to tissue dissociation, storage, and library preparation processes in scRNA-seq, during which immune cells may be preferentially preserved21. These observations suggest that the CosMx platform may offer an unbiased cellular map of liver tumors.

We distinguished malignant cells from non-malignant epithelial cells in the single-cell spatial data based on their transcriptomes and geographic locations in tumors (Supplementary Fig. 1E). Specifically, we performed clustering analysis of all epithelial cells and mapped each cluster to its spatial location. Clusters that mapped to tumor regions were labeled as malignant cells, whereas those mapped to adjacent non-tumor regions were defined as epithelial cells based on histological images (see Methods). We did not infer copy number variations from the transcriptome to confirm the malignancy of epithelial cells due to a limited number of targeted genes in the single-cell spatial data. In total, 792,192 malignant cells and 1,555,397 non-malignant cells were identified from the single-cell spatial data. Malignant and non-malignant cells from the scRNA-seq data have been determined from our previous analysis20.

Heterogeneous transcriptomic states in malignant cells
To understand the degree of ITH in PLC, we determined the transcriptomic states across all malignant cells from the single-cell spatial data using non-negative matrix factorization (NMF), a method has been successfully applied in defining tumor cell states based on scRNA-seq data4,5. We used scRNA-seq data as a reference to increase confidence in cellular states identification (Methods). Overall, we identified 12 distinct tumor cell states across all malignant cells tested, including 9 unique states of cell cycle, stress response, immune response and locomotion, metallothionein, EMT, hepatocyte-like, cholangiocyte-like, MHC-II, interferon response, as well as 3 mixed states of stress/metallothionein, stress/immune response and locomotion, and MHC-II/metallothionein (Fig. 2A, B, Supplementary Fig. 2A–C). These states were annotated based on the most highly expressed gene modules and their functional enrichment (Fig. 2B; Supplementary Data 3). Immune response and locomotion (CCL20, CXCL5, ICAM1) represented the most prevalent cellular states in malignant cells, highlighting the interactions between tumor cells and immune cells (Supplementary Fig. 2D). Other frequent states included stress/metallothionein, cholangiocyte-like, and cell cycle (Supplementary Fig. 2D). Interestingly, we observed a small proportion of MHCII-related tumor cells. Although MHCII molecules are primarily expressed by professional antigen-presenting cells, tumor cells can also express these molecules, and their specific expression by tumor cells has been linked to favorable outcomes in cancer patients22. In contrast to malignant cells, non-malignant epithelial cells were mainly enriched in metallothionein and stress-associated programs, which are related to the functions of hepatocytes (Supplementary Fig. 2D). We found most of the states were similar to those identified in pan-cancer studies, suggesting commonalities of tumor cell transcriptomic states among cancer types4,5. However, hepatocyte- and cholangiocyte-like states were uniquely observed in liver cancer, reflecting organ-specific cell lineages. Noticeably, each tumor contained a mixture of different tumor cell states, underscoring pronounced intratumor heterogeneity within each lesion (Supplementary Fig. 2E; Supplementary Data 4).

Landscape of non-malignant cells
We delineated the landscape of non-malignant cells by performing Harmony-based clustering of each cell type using the single-cell spatial transcriptome data23. In total, we identified 79 clusters, comprising 18 clusters of endothelial cells, 14 clusters of fibroblasts, 16 clusters of myeloid cells, 10 clusters of B cells, and 21 clusters of T cells, reflecting a complex and heterogeneous landscape of immune and stromal cells in liver cancer (Fig. 2C and Supplementary Fig. 3). We further annotated these clusters into 44 cellular states based on differentially expressed genes (Fig. 2C and Supplementary Fig. 3). For example, in the fibroblast population, matrix fibroblasts (COL1A1), antigen-presenting fibroblasts (CD74, IGKC), smooth muscle cells (RGS5), EMT-like fibroblasts (KRT17), tumor-like fibroblassts (VEGFA), CXCL12+ CAFs (CXCL12, IGF2), IL32+ CAFs (IL32, CXCL10), and CCL19+ CAFs (CCL19) were determined. Similarly, different cellular states in endothelial cells, myeloid cells, and lymphoid cells were defined (Fig. 2C). We validated cluster robustness by comparing them with clusters derived from individual patients and found that most patient-derived clusters mapped predominantly to specific Harmony-based clusters, indicating minimal over-correction by Harmony (Supplementary Fig. 4). To determine if the non-malignant cell subtypes resolved from the single-cell spatial data and the scRNA-seq data were consistent, we compared the clusters in each major cell type derived from the two approaches. We found that most clusters from the spatial data could be well-matched to specific clusters in the scRNA-seq data, indicating that the spatial approach effectively captures major cellular information despite targeting fewer genes than scRNA-seq (Supplementary Fig. 5). Interestingly, we observed some unique clusters in the spatial data. For example, in myeloid cells, the TUBB+ proliferative tumor-associated macrophages (TAMs) were only found in the spatial data, suggesting that the single-cell spatial method may identify potential novel cell clusters (Supplementary Fig. 5)24,25. Collectively, the single-cell spatial transcriptomic profiles resolved by CosMx offer a detailed landscape of non-malignant cells in liver cancer.

Spatial distribution of cells within and across tumor regions
To assess intratumor heterogeneity, we analyzed the spatial distribution of cell (sub)types within and across tumor regions. Overall, epithelial cells and T cells were markedly decreased in both the tumor border and core relative to adjacent non-tumor tissue, whereas myeloid cells and fibroblasts gradually increased from non-tumor region to the border and core (Supplementary Fig. 6A). Among epithelial cells, malignant cells increased and non-malignant epithelial cells decreased from the non-tumor region toward the tumor core (Supplementary Fig. 6A). Detailed comparison of tumor cell states revealed largely comparable proportions between the border and core regions in our cohort (Supplementary Fig. 6B). While most non-malignant cell subtypes showed no significant differences across these regions, several displayed region-specific enrichment. At the tumor border, two immune subtypes (Kupffer cells and CD8+ effector T cells) and four stromal subtypes (liver sinusoidal endothelial cells [LSECs], post-capillary venule endothelial cells, CCL19+ CAFs, and CXCL12+ CAFs) were significantly enriched compared to the tumor core (Supplementary Fig. 6C). This observation was consistent with the typical localization of Kupffer cells and LSECs in normal liver tissue. In contrast, tip cells and smooth muscle cells were elevated in the tumor core, along with reg-TAMs, suggesting a potential pro-angiogenic and immunosuppressive microenvironment (Supplementary Fig. 6C)26–29. We also investigated the spatial variation of cells within each tumor region by calculating Moran’s I to quantify spatial autocorrelation. For all major cell types, Moran’s I values were positive, indicating that each cell type exhibited some degree of spatial autocorrelation within individual tumor regions (Supplementary Fig. 7A). Notably, Moran’s I increased for epithelial cells, myeloid cells, fibroblasts, and endothelial cells in the tumor border and core compared with adjacent non-tumor regions, suggesting stronger spatial clustering of these populations within tumors (Supplementary Fig. 7A). In contrast, Moran’s I values for T cells and B cells decreased in the tumor border and core relative to non-tumor regions, reflecting a more dispersed distribution within the tumor (Supplementary Fig. 7A). These observations are consistent with the clustered distribution of T cells and B cells around portal veins in normal liver, whereas other cell types are more evenly dispersed in normal tissue30. Malignant cell states displayed low but positive Moran’s I values, with no significant differences between border and core, suggesting only marginal spatial autocorrelation of each malignant cell state within individual tumor regions for the samples in our cohort (Supplementary Fig. 7B).

Characterization of the surroundings of individual malignant cells by defining spatial dynamics network
A successfully evolved tumor represents a spatially well-organized ecosystem, where malignant cells actively interact with their local environment. This dynamic interplay continuously shapes tumor cell functions, fueling ITH. To uncover the roles of spatial context in driving ITH, it’s essential to delineate the surroundings of each individual malignant cell. The resolved single-cell spatial data in this study provide a unique and powerful approach for this purpose. We characterized the surroundings of each malignant cell by determining the composition of its neighbors within 40 µm, as in our previous study (Fig. 3A)31. We further performed clustering analysis of all malignant cells based on their surroundings using the Louvain algorithm (Fig. 3A). Consequently, malignant cells falling into the same cluster shared similar spatial neighbors. This strategy is independent of transcriptome-based clustering, where tumor cells within a certain cluster have similar transcriptomes. We tested the stability of the malignant cell clusters by applying different radii of 20 µm, 60 µm, 80 µm, and 100 µm, where similar results were found to those derived from 40 µm, demonstrating stable clusters of malignant cells using the neighborhood information (Supplementary Fig. 8A, B). In total, 14 clusters of malignant cells (corresponding to 14 types of surrounding scenarios) with distinct surrounding environments were determined (Fig. 3B, C, Supplementary Fig. 8C, D). We annotated different types of surrounding environments based on their compositions and designated each type as a Spatial Dynamics Network (SDN) around malignant cells (Fig. 3A, C). Vascularized tumors represented a major group of SDNs (SDN-T1–SDN-T5), with endothelial cells in the surroundings of malignant cells. We also identified a tumor-dominant SDN (SDN-T8), in which malignant cells comprised the majority of surrounding cells. Additionally, we found five SDNs (SDN-T10–SDN-T14) with an enrichment of myeloid cells and fibroblasts, as exemplified in an iCCA sample (Fig. 3C, D). These findings highlight the diverse microenvironments surrounding malignant cells in liver cancer.

Tumor transcriptomic states are related to their spatial contexts
To understand the roles of spatial context in shaping ITH, we analyzed the relationship between tumor cell transcriptomic states and the SDNs. Specifically, we performed an enrichment analysis of malignant cell states on the SDNs using a random shuffling strategy (Method). Notably, tumor cell states displayed non-random spatial distributions, with each state preferentially enriched in specific SDNs (Fig. 4A). For example, EMT-like tumor cells were enriched in SDNs (SDN-T1–SDN-T5) that were related to vascularized tumors. Consistently, an increased density of endothelial cells along with EMT transition has been reported32. In addition, we found EMT-like tumor cells were in significantly closer spatial proximity to endothelial cells compared to the tumor cell states that were enriched in completely different SDNs (Fig. 4B). We observed unique populations of IGF2+ lipid-associated macrophages (LA-TAM) in SDN-T1 and SDN-T2, where EMT was the most strongly enriched (Figs. 3C and 4A). IGF2 has been demonstrated to be secreted by macrophages to regulate EMT33. We also found CXCL12+ CAFs expressing IGF2 and CXCL12 in SDN-T1 and SDN-T2, suggesting a pivotal role of the local environment in regulating EMT phenotype of tumor cells (Fig. 3C). In contrast to the complex local microenvironment of EMT-like tumor cells, cell cycle-related tumor cells were strongly enriched in tumor-dominant SDNs (SDN-T7–SDN-T9) with few stromal and immune cells. This was further supported by cell density analysis, with a higher density of tumor cells surrounding cell cycle-related tumor cells than others (Fig. 4C). Additionally, we noticed that cell cycle- and EMT-related tumor cells were enriched in completely different SDNs. It’s well known that tumor cells going through EMT have a decreased ability to proliferate34. We further demonstrated distinct local environments of the two tumor cell states. Interferons are cytokines that help trigger the immune system to eradicate pathogens or tumor cells35. We observed an enrichment of interferon response-related tumor cells in SDN-T10 and SDN-T2, where myeloid cells and CD8+ T cells were found. These tumor cells expressed markers of IFIT1 and ISG20, which are related to antiviral defenses36,37. Noticeably, interferon response-related tumor cells were mainly found in a patient with infections of the hepatitis C virus and human immunodeficiency virus in our cohort, suggesting a potential antiviral response in this patient (Supplementary Fig. 2E; Supplementary Data 1). We also found a strong enrichment of cholangiocyte-like tumor cells and immune response and locomotion-related tumor cells in the SDNs with fibroblasts and myeloid cells (SDN-T11–SDN-T13). Ligand-receptor interaction analysis between tumor cell states and their local environments demonstrated enriched communication patterns in each state, further supporting the spatial preference of different tumor cell states (Supplementary Fig. 9). For example, interactions mediated by VEGF were mainly observed between EMT-related tumor cells and their surrounding environments (Supplementary Fig. 9).

To validate the spatial preference of tumor cell states, we analyzed liver tumor samples from three publicly available datasets based on the 10× Visium platform6,38,39. Since 10× Visium data is at a spot (55 µm in diameter) resolution, the methods designed for single-cell spatial data cannot be directly applied to determine cellular states and SDNs here. Therefore, we developed pseudo-bulk gene signatures for each tumor cell transcriptomic state and SDN by averaging the associated genes from our single-cell spatial data (Methods). We further applied these gene signatures to the three datasets and used a correlation-based approach to determine the enrichment of each tumor cell state in the SDNs (Methods). We observed highly consistent enrichment patterns with those found in the single-cell spatial data, with a much stronger enrichment of tumor cell states in SDNs that showed significant association in the single-cell spatial analysis compared to other SDNs (Fig. 4D). Collectively, these observations suggest that the transcriptomic states of tumor cells are not randomly distributed within a tumor. Instead, they are closed associated with their surrounding microenvironments, underlying the potential roles of spatial context in driving ITH.

Tumor cell villages are identified using graph attention networks
Diverse tumor cells may spatially coordinate to collectively drive tumor growth and progression, conferring survival advantages to the tumor and enhancing its resilience against interventions40. We define this spatial coordination of diverse tumor cell transcriptomic states, supported by specific local environments, as tumor cell “villages.” Like human villages, tumor cell “villages” may use communal mechanisms to promote growth and enhance defense, reducing individual cell vulnerability.

To determine tumor cell villages, we used graph attention networks, which are deep learning neural networks specifically designed for analyzing graph-structured data41. These networks leverage an attention mechanism to enhance feature learning by weighing the importance of neighboring nodes, allowing us to capture the complex spatial relationships between diverse tumor cell states. In general, malignant cells within a tumor were represented as connected graphs, where nodes indicated individual cells characterized by their cellular states and SDNs, and edges linked nodes within a 40 µm distance, consistent with the SDN analysis (Fig. 5A). The embeddings derived from the graph attention networks were then applied for clustering analysis to define tumor cell villages (Fig. 5B). We excluded three small clusters with fewer than 5000 cells each, accounting for less than 1% of total malignant cells. Eight distinct tumor cell villages were then identified (Fig. 5C and Supplementary Fig. 10A). Village 1 and Village 2 contained a mix of tumor cell states with a vascularized environment. By contrast, Village 4 and Village 7 were dominated by cell cycle-related tumor cells, akin to cancer germinal centers. Villages 3, 5, and 6 were enriched in cholangiocyte-like and immune response and locomotion-related tumor cells, surrounded by a microenvironment of fibroblasts and macrophages, reflecting a well-developed ecosystem of tumor cells with ample resources. Village 8 was primarily associated with immune response and locomotion. As an example, 1CB (an iCCA sample) featured Villages 3, 5, and 6 in distinct geographic locations (Fig. 5D). We observed that tumor cell village compositions were more consistent across regions within the same patient than across different patients in our cohort (Supplementary Fig. 10B).

We determined surrogate markers for tumor cell villages through differential gene expression analysis (Supplementary Fig. 10C; Supplementary Data 5). The identified genes were found largely specific to each tumor cell village and can accurately predict the presence of their corresponding villages (Supplementary Fig. 10D–F), suggesting that these genes may effectively serve as surrogates for defining tumor cell villages. Using these genes, we validated the identified tumor cell villages with spatial transcriptome data from 46 HCC or iCCA tumors profiled on the 10X Visium platform6,38,39,42 (Supplementary Fig. 11A, B). To examine how tumor cell villages relate to patient prognosis, we performed survival analysis on 37 HCC or iCCA patients with single-cell transcriptome data1, stratifying the patients by surrogate gene-based hierarchical clustering. We found that patients enriched for Villages 3–8 genes signatures had significantly poorer overall survival compared to those with features related to Villages 1 and 2 (Fig. 5E and Supplementary Fig. 11C). A consistent trend was observed when restricting the analysis to HCC patients alone (Supplementary Fig. 11D). We didn’t analyze iCCA patients separately, as all were predominantly enriched in Villages 3–8. The survival associations were further confirmed using bulk transcriptomic data from 674 HCC patients in the TCGA, LCI (Liver Cancer Institute)43, and Mongolia44 cohorts, where consistent links between tumor cell villages and clinical outcomes were observed (Fig. 5F–H and Supplementary Fig. 11E). These analyses highlight the presence of well-coordinated tumor cell villages in liver cancer and suggest that tumor cell villages are linked to patient outcomes.

To investigate the mechanisms underlying the observed clinical associations, we performed gene set enrichment analysis. Villages 1 and 2 were enriched in metabolic-related pathways, such as bile acid metabolism, adipogenesis, and fatty acid metabolism (Supplementary Fig. 11F). In contrast, pathways such as MYC targets, E2F targets, G2/M checkpoint, mitotic spindle, reactive oxygen species pathway, PI3K/AKT/mTOR signaling, and angiogenesis were observed across multiple Villages within 3–8, corresponding to programs that promote cell growth, survival, resistance to cell death, and neovascularization (Supplementary Fig. 11F). Although a small subset of EMT-like cells was detected in Village 2, EMT pathways were not enriched. Notably, while Villages 1 and 2 contained more endothelial cells, the endothelial cells in Villages 3–8 exhibited significantly higher expression of tumor-specific endothelial signatures, consistent with the angiogenesis pathway observed in several of these villages (Supplementary Fig. 11G)45. We did not observe significant differences in the overall abundance of inflammatory-related fibroblasts (CXCL12+ CAFs, IL32+ CAFs, and CCL19+ CAFs) between Villages 1–2 and Villages 3–8 (Supplementary Fig. 11H). Collectively, these results suggest that Villages 3–8 may represent more aggressive tumor ecosystems compared with Villages 1 and 2, which may underlie their unfavorable clinical outcomes.

Molecular co-dependencies linked to tumor cell villages
To uncover key factors driving cell coordination within each tumor cell village, we developed a method to evaluate the molecular co-dependencies between tumor cells and non-tumor (stromal or immune) cells. Specifically, we calculated the correlations between each pair of genes from tumor and non-tumor cells based on a 40 µm distance, with an assumption that highly correlated gene pairs are strongly co-dependent (Methods). Using this method, we identified the top correlated gene pairs from each village and found that these gene pairs are village-specific, with significantly higher correlation scores in their respective villages compared to others (Fig. 6A; Supplementary Data 6). We further evaluated the gene-pair correlations in each village by increasing the spatial distance thresholds from 40 µm to 80 µm, 160 µm, 320 µm, and by extending beyond the village boundaries. Notably, a gradual decline in correlation scores was observed as the spatial separation of cells increased, underscoring the spatial specificity of the gene-pair relationships (Fig. 6B). To demonstrate the importance of these gene-pair relationships to each tumor cell village, we employed a random forest model to predict village identities of the tumor cells based on the identified gene pairs. We observed a significant decline in prediction accuracy when the non-tumor cells surrounding the tumor cells within each village were randomly shuffled, compared to the accuracy achieved with the original village configurations (Fig. 6C). These findings suggest that molecular co-dependencies may be crucial for maintaining the integrity of each village, and their perturbation may potentially destabilize or even collapse the village structure.

Among the top-ranked molecules from tumor cells, we focused on serine protease inhibitor Kazal-type 1 (SPINK1) due to its unique expression in tumor cells (Fig. 6D and Supplementary Fig. 12A). SPINK1-related gene pairs were predominantly found in Villages 6, with a few detectable in Village 3 (Fig. 6E). Most of the SPINK1-associated genes were matrix CAFs (mCAFs)-related genes, including COL1A1, COL1A2, COL5A1, COL4A1, BGN, and MGP (Fig. 6F; Supplementary Data 6). This spatial relationship was validated using Visium transcriptome data, which showed that SPINK1+ spots were significantly closer to the paired gene-enriched spots than expected by random chance (Fig. 6G). Multiplex immunofluorescence for SPINK1 and COL1A1 further confirmed the colocalization of the two proteins in patient 3C, who had the highest proportion of Village 6 among all patients (Supplementary Figs. 10B and 12B). Gene set enrichment analysis suggested that SPINK1+ tumor cells were associated with metastasis-related features, while SPINK1- tumor cells were enriched in metabolism- and coagulation-related pathways (Fig. 6H). Additionally, single-cell transcriptome analysis demonstrated that SPINK1+ tumor cells engaged more extensively with mCAFs than other CAFs through ligand-receptor interactions (Fig. 6I; Supplementary Data 7 and 8). SPINK1 itself was identified as one ligand in these interactions with EGFR serving as its receptor (Supplementary Data 7)46,47. Consistently, SPINK1 downstream genes were expressed at significantly higher levels in mCAFs compared with other CAFs (Supplementary Fig. 12C). Similarly, the EGFR downstream genes were significantly upregulated in mCAFs (Supplementary Fig. 12D). When the relationship of SPINK1 and its paired genes were perturbed by randomly shuffling non-tumor cells, we found that the prediction accuracy of village identities by the random forest model declined significantly compared with the original village configurations, with a greater decrease observed in Village 6 than in Village 3, consistent with the higher enrichment of SPINK1-related gene pairs in Village 6 (Fig. 6J). The prediction accuracy was further declined when all the identified correlated genes in these villages were perturbed (Fig. 6J). These observations suggest the molecular co-dependencies between SPINK1+ tumor cells and mCAFs may play an important role in sustaining the integrity of Villages 3 and 6. Collectively, these findings indicate that the defined tumor cell villages may reveal unique molecular co-dependencies among cells within each village, opens avenues for developing treatment strategies aimed at disrupting these co-dependencies to destabilize the tumor.

Methods
Humans sample collection
This study included seven liver cancer patients (four HCC, three iCCA; four males, three females) from the University Medical Center in Mainz and the NIH Clinical Center. Detailed information of this cohort can be found in Supplementary Data 1. The ages of the patients ranged from 63 to 77, with a median age of 72. A total of 50 samples were collected from these patients. We conducted single-cell spatial transcriptome profiling of 15 samples. Additionally, we generated single-cell transcriptome data for the remaining samples20. Sample collection was performed with informed consent from patients. This study was approved by the ethics committee of the University Medical Center in Mainz and the National Institutes of Health.

CosMx SMI sample preparation
Formalin-fixed, paraffin-embedded (FFPE) tissue sections were prepared for CosMx SMI profiling following the procedures described in He et al.19. Briefly, five-micron tissue sections were cut and placed on Leica BOND PLUS slides. Slides were baked overnight at 60°C to improve tissue adherence to the slide. Tissues were deparaffinized in Xylene twice for 5 min each. Then, tissues were successively dehydrated in 100% ethanol twice for 5 min each and subsequently dried at 60°C for 5 more min. Tissues were then prepared for in-situ hybridization (ISH) by heat-induced epitope retrieval (HIER) at 100°C for 15 min using 1× Target retrieval solution (NanoString CMx Slide Prep Kit, FFPE RNA 121500006) in a Bio SB cooker. After target retrieval, tissue sections were rinsed with diethyl pyrocarbonate (DEPC)-treated water (DEPC H2O, ThermoFisher), washed in ethanol for 3 min and dried at room temperature for 30 min. Once the slides dried, an incubation frame was placed over the tissue (CosMx FFPE Slide Prep Kit RNA). Tissue was then digested with Proteinase K (provided by NanoString) 3 µg/ml at 40°C for 30 min. Tissue sections were rinsed twice with DEPC H2O, incubated in 1:400 diluted fiducials (Bangs Laboratories) in 2× saline sodium citrate and Tween (0.001% Tween 20, Teknova) for 5 min at room temperature, and washed with 1× PBS (ThermoFisher) for 5 min. After digestion and fiducial placement, tissue sections were fixed with 10% neutral buffered formalin (NBF) for 1 min at room temperature. Fixed samples were rinsed twice with Tris-glycine buffer (0.1 M glycine, 0.1 M Tris-base in DEPC H2O) and once with 1× PBS for 5 min each before blocking with 100 mM N-succinimidyl (acetylthio) acetate (NHS-acetate, ThermoFisher) in NHS-acetate buffer (0.1 M NaP, 0.1% Tween PH 8 in DEPC H2O) for 15 min at room temperature. The sections were then rinsed with 2× saline sodium citrate for 5 min. NanoString ISH probes were prepared by incubation at 95°C for 2 min and placed on ice, and the ISH probe mix (1 nM 980 plex ISH probe, 10 nM Attenuation probes, 1× Buffer R, 0.1 U/µL SUPERase•In™ [Thermofisher] in DEPC H2O) was pipetted on the tissue slides. An adhesive SecureSeal Hybridization Chamber (Grace Bio-Labs) was placed over the tissue. The hybridization chamber was sealed to prevent evaporation, and hybridization was performed at 37°C for 16–18 h. Tissue sections were rinsed of excess probes in 2× SCT for 1 min and washed twice in 50% formamide (VWR) in 2× SSC at 37°C for 25 min, then twice with 2× SSC for 2 min at room temperature. Samples were then prepared for antibody cocktail incubation. First, slides were incubated with nuclear stain 1:40 (provided by NanoString) for 15 min at room temperature, protected from light. After a short PBS wash, tissues were incubated with a 4-fluorophore-conjugated antibody cocktail against CD298/B2M (488 nm), PanCK (532 nm), CD45 (594 nm), and CD68 (647 nm) proteins and DAPI stain in the CosMx SMI instrument. Antibodies were also included in the CosMx slide preparation kit by NanoString. After 1 h incubation at room temperature protected from light, slides were washed three times with 1× PBS 5 min each. A custom-made flow cell was affixed to the slide in preparation for loading onto the CosMx SMI instrument. Flow cells were then stored in SSC while the instrument was being prepared for loading the samples.

CosMx imaging
RNA target readout on the CosMx SMI instrument was performed according to the procedures described by He et al.19 Briefly, the assembled flow cell was loaded onto the instrument, and Reporter Wash Buffer was flowed to remove air bubbles. A preview scan of the entire flow cell was taken, and fields of view (FOVs) were placed on the tissue to match regions of interest identified. In our case, we selected as many FOVs as needed to completely cover the whole tissue on the slide. RNA readout began by flowing 100 µl of Reporter Pool 1 into the flow cell and incubation for 15 min. Reporter Wash Buffer (1 mL) was flowed to wash unbound reporter probes, and Imaging Buffer was added to the flow cell for imaging. Eight Z-stack images (0.8 µm step size) for each FOV were acquired, and photocleavable linkers on the fluorophores of the reporter probes were released by UV illumination and washed with Strip Wash buffer. The fluidic and imaging procedure was repeated for the 16 reporter pools, and the 16 rounds of reporter hybridization-imaging were repeated multiple times to increase RNA detection sensitivity. Then, tissue and cell membrane morphology visualization was performed using oligonucleotide-conjugated antibodies hybridized to a specific pool of barcode reporters. After unbound antibodies and DAPI stain were washed with Reporter Wash Buffer, Imaging Buffer was added to the flow cell, and eight Z-stack images for the 5 channels (4 antibodies and DAPI) were captured.

CosMx SMI data quality-control and preprocessing
With the coordinates and the raw transcript counts of individual cells, we filtered cells based on the following criteria: ≥20 counts detected and ≥10 genes detected; <10% negative probes and ≤5 negative probe counts. Cells with outlier cell area (p-value < 0.01) based on Grubb’s test were also removed for downstream analysis. A total of 2,347,589 cells passed the quality control with a median of 93 features and 174 counts per cell. Data normalization was performed using the R package Seurat (v5.1.0).

Major cell type annotation
Major cell type annotation (T cells, B cells, myeloid cells, endothelial cells, fibroblasts, and epithelial cells [including hepatocytes, cholangiocytes, and tumor cells]) of the single-cell spatial data was based on label transfer using the scRNA-seq dataset20. The R package Harmony (v1.2.0) was used to integrate the CosMx SMI data with the scRNA-seq data. The top 20 principal components (PCs) were used for the integration, and the top 20 harmony embeddings were used to cluster the cells in the integrated dataset with default parameters. For each derived cluster, the percentage of different cell types from the scRNA-seq dataset was calculated. Clusters were annotated based on the highest proportion of cell types from scRNA-seq. A cluster was considered as unclassified if the highest proportion < 50% or the ratio of the second-highest proportion against the highest proportion > 50%. The analyses were performed for each sample separately. The expression score for each cell type was calculated as the mean values of the cell type-specific genes (Supplementary Data 2).

Malignant cells were determined based on their transcriptomes and geographic locations in a tumor. For epithelial cells from each individual patient, we identified different clusters of epithelial cells using Seurat. Specifically, after data normalization, most variable features were selected based on mean expressions and standard variances of genes. Top 20 PCA components were used for detecting neighbors, and Louvain clusters of epithelial cells were determined using default parameters. Each cluster was subsequently mapped to its spatial location. Clusters of epithelial cells located within tumor regions, based on the histological images, were classified as malignant cells, while those in non-tumor regions were annotated as non-tumor epithelial cells.

Consistency of CosMx SMI and scRNA-seq data
To compare the gene expression profiles between CosMx SMI and scRNA-seq data generated from the same patients, we grouped the cells into six major cell types, including T cells, B cells, myeloid cells, endothelial cells, fibroblasts, and epithelial cells. We selected the genes that were detected in both datasets and calculated the average expression of each gene in individual cell types. We further calculated the Pearson correlation coefficient (R function “cor.test”) between averaged gene expression profiles from CosMx SMI data and scRNA-seq data to measure the consistency of the two profiling approaches.

Define tumor cell transcriptomic states based on gene modules
Due to the limited number of genes in the single-cell spatial data compared with scRNA-seq data, we determined gene modules in malignant cells based on scRNA-seq data and further applied those identified modules into the single-cell spatial data to define tumor cell transcriptomic states. Gene module identification was conducted based on a previous published non-negative matrix factorization (NMF)-based method5. We included two liver cancer single-cell datasets to improve the robustness and sensitivity of gene module identification15,20. For each sample, the top 500 highly variable genes (HVGs) were used for data scaling and the NMF analysis. The number of cells was down-sampled to 1000 for the samples containing over 1000 tumor cells. Negative values were set to 0, and genes with expression values of 0 across all cells within each sample were discarded. We then performed NMF on the scaled gene expression matrix for each sample separately using the R package NMF (v0.27). The “nsNMF” method and “nndsvd” seeding method were applied. The rank was set to 10 if the number of malignant cells per sample was less than 400, and to 30 for samples containing 400 or more malignant cells. The contribution matrices of genes in different modules for each sample were obtained. A gene was assigned to a module if it showed the highest contribution to that module across all modules within a sample. Modules obtained from individual samples were further filtered by the Jaccard index. Gene modules with ≥5% overlap with at least two other modules were retained. For each paired gene, the number of individual tumor modules in which they co-existed (designated as the connection value) was calculated, and a gene–gene connection matrix was constructed for gene pairs. Genes with a connection value greater than three with at least four other genes were retained for downstream analysis. The graph was clustered using the infomap clustering method from the R package igraph (v2.0.3). The final modules with at least 5 genes were retained. Using this approach, we identified 11 modules across all samples. We annotated the gene modules based on gene ontology (GO) terms from MSigDB55, tumor-derived signatures from Neftel et al.56, Puram et al.57, Ji et al.58, and Barkley et al.5, Hypergeometric distribution test was used to determine the enrichment of our identified consensus module genes in annotated gene sets with R function “phyper”. FDR method was used to adjust the p-values.

Based on the gene modules identified from scRNA-seq data, module scores of individual epithelial cells from CosMx SMI dataset were calculated using the method from Barkley et al.5, For each module, 10,000 randomized gene lists were generated, each containing the same number of genes as the module. The average gene expression for each random list and the module was then calculated. The score for the module was defined as: score = log10(n+1)/10000, where n represents the number of random gene sets with higher average expression than the module. The score was further linearly rescaled to between 0 and 1. We excluded “Metabolism” and “Apoptosis” modules from further analysis due to less than two genes detected in the CosMx SMI dataset. We removed cells with a score variance less than 0.15 and a mean score less than 0.2 across all modules to increase the confidence of module determination. The derived matrix was used to perform PCA and UMAP analyses. The epithelial cell clusters were determined by a graph-based method using the top 8 PCA components with default parameters in Seurat (v5.1.0). Average gene module scores were calculated for each cluster. Clusters with similar module scores were merged based on their hierarchical relationship, and a total of 32 epithelial clusters were obtained. Epithelial clusters were annotated to specific transcriptomic states based on the highest gene module scores.

Cell subtypes in non-malignant cells
For each non-malignant cell type (myeloid cells, T cells, B cells, endothelial cells, and fibroblasts) in the single-cell spatial transcriptome data, Harmony (v1.2.0) was used to integrate cells from different patients. Top 20 PCA components were used to perform harmony analysis, and the top 20 harmony embeddings were used for UMAP and clustering. Clusters were annotated based on differentially expressed genes with thresholds of adjusted p-value < 0.05 and log2 fold change ≥ 0.5. Marker genes from literature were used as reference in this analysis, including markers from Goveia et al.27 for endothelial subtypes, from Ma et al.29 for myeloid subtypes, from Cords et al.59 and Zhang et al.60 for fibroblast subtypes, from Morgan et al.61, Yang et al.62, and Zhang et al.63 for B cell subtypes, and from Guo et al.64, Zhang et al.65, Zheng et al.66, and Ma et al.1 for T cell subtypes.

Align cell subtypes between CosMx SMI data and scRNA-seq data
To compare cell subtypes identified from the CosMx SMI data and the scRNA-seq data, we performed correlation analysis of the clusters identified from the two datasets. Specifically, we conducted cluster analysis of each non-epithelial cell type, i.e., myeloid cells, T cells, B cells, endothelial cells, and fibroblasts, from the scRNA-seq data with top 10 PCA and default resolutions. We selected the common genes from the CosMx SMI and scRNA-seq data to calculate the Pearson correlation coefficients between each pair of clusters from the two datasets. A CosMx™ SMI cluster was considered present in the scRNA-seq data if its correlation coefficient with any scRNA-seq cluster exceeded 0.5.

Characterization of the surroundings of individual malignant cells
To characterize the surroundings of individual malignant cells, we determined the composition of surrounding cells for each targeted malignant cell in the CosMx™ SMI data. The surrounding cells are defined as those located within a 40 µm radius of the targeted cell. We included three layers of information to characterize the surrounding cell compositions, including major cell types, cell states/subtypes, and cell clusters. This vector of cell compositions was used as features to describe the targeted cell. We further performed PCA analysis of the surrounding composition matrix for all malignant cells and used the top 10 PCA components to generate a UMAP representation with default parameters. Clusters of malignant cells were determined using the Louvain algorithm at a resolution of 0.6. Each cluster of malignant cells was annotated based on the composition of the surrounding cells.

To determine the distance for surrounding cell characterization, we tested the radii of 20, 40, 60, 80, and 100 µm. We performed pairwise comparison of the derived results using a normalized mutual information score. We also tested the significance of the comparisons based on a random shuffling strategy. We used 1000 times of randomly shuffled malignant cells as a background to calculate the p-value. Compared with 20 µm, malignant cell clusters are more stable with a larger radius. At a radius of 40 µm, the median number of surrounding cells for all malignant cells is 33. Considering that a larger radius will introduce a higher number of surrounding cells with a higher probability of noise to characterize a targeted cell, we selected a radius of 40 µm to determine the surroundings of individual malignant cells in our liver cancer single-cell spatial data.

Enrichment of tumor cell states on SDNs
In the CosMx™ SMI data, the enrichment of tumor cell states on SDNs was calculated as the proportions of each tumor cell state with a specific SDN. We used a random shuffling approach (1000 times of randomly shuffling of tumor cell states) to calculate p-values and to determine the significance of the enrichment. To identify significantly enriched tumor cell states, p-values were calculated in a stringent way as p-value = n/1000, where n represents the number of cases with expected scores ≥ 0.8*observed score. A p-value less than 0.05 was considered significant.

We used 39 samples from three publicly available 10× genomics Visium datasets from, e.g., Liu et al.38, Wu et al.6, and Zhang et al.39 to validate the enrichment of tumor cell states on SDNs. Among them, 24 samples from the tumor core were included in our analysis. To determine tumor cell state for each spot, we generated reference pseudo-bulk gene signatures based on the CosMx™ SMI data. We averaged gene expression of the same tumor cell states to build an expression vector for each state. Genes with mean expression < 0.1 or standard deviation < 0.1 across all cell states were removed to retain the variable genes. We further selected the genes that were detected in the three datasets to ensure consistency. For each Visium dataset, Spearman correlations were then calculated between reference cell state signatures and the gene expression of each spot from the Visium data, resulting in a correlation matrix (Matrix I) of cell states and spots. High correlation values indicate a high probability of the spots containing tumor cells with the indicated cell states. To determine the SDNs for each spot in the 10× Visium datasets, we collected all cells related to each SDN in the CosMx data. The averaged gene expression related to each SDN was determined. Similar to the analysis of tumor cell states, we used the same criteria to filter genes and calculate correlations to determine the SDN of each spot. A correlation matrix (Matrix II) of SDNs and spots was generated. To infer the tumor cell state enrichment on SDNs, we calculated the Pearson correlation between tumor cell states (from Matrix I) and SDNs (from Matrix II) across all spots. The Pearson correlation coefficients were linearly rescaled to a range between 0 and 1.

Physical distance between malignant cells to endothelial cells
For each malignant cell, the distance to the closest endothelial cell was calculated using the R function “kNN” from dbscan (version 1.1-12). R function “kdtree” was used to search the nearest cell.

Identification of tumor cell villages
To identify tumor cell villages, we constructed an undirected graph for tumor cells from each sample, where nodes represent tumor cells and edges indicate connections of tumor cells within 40 µm distance. The cell state and the SDN were represented by one-hot encoding and combined as node features. Graph Attention Network (GAT) model implemented in PyTorch (v1.21.1) was applied to train the connected graphs for node clustering. In contrast to other methods, GAT utilizes attention mechanisms. During training, the model learns the importance of each neighboring node’s features to the target node as attention coefficients for the target node. The coefficients are then normalized by the SoftMax function to sum to 1. We used this attention scoring as a representation of the interaction strength between nodes. These coefficients are used during feature aggregation in that the learned node feature is a weighted sum of neighbor information according to the attention coefficients. A 128-dimensional hidden representation was used for message passing between neighbors. Three convolutional layers was used in our analyses to aggregate neighborhood information of 3-hops. The Exponential Linear Unit (ELU) activation function was applied for each layer except for the last layer, followed by a 25% dropout rate for regularization. The model outputted a 15-dimension feature vector for each cell. To maximize the similarity of the embedded neighboring cells while minimizing that of distant cells, we used an evaluation method by Shiao et al.67 Specifically, a mini-batch of 10,000 target cells was randomly sampled from the input graph. For each target cell, we sampled two types of cells: one directly connected to the target cell, and another chosen from a random position in the graph, with the latter assumed to be spatially unrelated to the target cell. The evaluation function maximizes similarity between neighboring cells and minimizes similarity of unrelated cells. The Adam optimization algorithm was used during training with a learning rate of 1e-5 for 100 epochs. We performed dimensional reduction on the derived embeddings from GAT for all the tumor cells using PCA. The top 10 PCA components were used for Louvain clustering at a resolution of 0.5, resulting in 10 tumor clusters. Clusters 9 and 10 were discarded from the downstream analysis due to a small number of cells, accounting for less than 1% of total malignant cells.

Village graph of malignant cells
Village graph represents the compositions of tumor cell states and SDNs associated with each village. The sizes of the dots indicate the proportion of tumor cell states in each village. Tumor cell states with a proportion of at least 5% were included in the graph. Interaction weights between tumor cells were obtained from the graph attention auto-encoder model described above, with higher values indicating stronger interactions. The interactions were visualized using the Fruchterman-Reingold layout. The outer ring stands for the SDNs, with different areas indicating the proportions of different SDNs related to each village.

Define tumor cell villages using gene signatures
We determined gene signatures for each village based on differential gene expression analysis of the CosMx™ SMI data. Genes with fold changes ≥ 1 and p-values < 0.05 were selected as surrogate signatures for each village. To determine the specificity of the signatures to each village, we constructed a pseudo-bulk dataset based on the CosMx™ SMI data. For each tumor cell, surrounding cells within a 40 µm radius, measured by Euclidean distance, were collected. Average gene expression across cells was then calculated to generate pseudo-spot bulk gene expression. For each pseudo-spot, we determined the village scores as the mean expression of each village-related genes. Each pseudo-spot was assigned to the village with the highest village score. A confusion matrix was calculated between the predicted villages and ground truth from CosMx™ SMI data analysis. Accuracy, precision, recall, and F1 score were calculated to measure the performance of our village detection model. F1 score was calculated as follows: F1 = 2 * (precision * recall) / (precision + recall).

For the 10× Visium dataset, we used the same method to determine tumor cell villages for the spots. For each dataset, village scores were then linearly rescaled to a range between 0 and 1 across all spots. The confidence of village assignment for each spot was calculated as the difference between the highest village score and the second-highest village score. For the scRNA-seq data, pseudo-bulk gene expression for each sample was determined, followed by tumor cell village score calculation.

Overall survival
For each cohort, tumor samples were clustered based on the village scores. Hierarchical relationship of the samples was constructed using Pearson correlation as distance and “complete” agglomeration method. Samples in each cohort were grouped into two clusters based on hierarchical relationships. Kaplan-Meier plots were generated and visualized using “Surv” function from R package Survival (v3.7-0).

Cell-cell interactions based on single-cell spatial transcriptome data
To infer cell-cell communications from the CosMx™ SMI data, we used the COMMOT method68, with CellChatDB as the ligand-receptor (LR) database. For each sample, we included only LR pairs identified in ≥ 100 cells or in ≥ 5% of all tested cells. Cells within a 50 µm distance were considered for interaction calculations. To identify LR interactions associated with tumor cell states, we used a multi-step approach. First, LR scores were calculated for each tumor cell state using tl.cluster_communication function in COMMOT, with cells categorized into B cells, T cells, endothelial cells, epithelial cells, fibroblasts, myeloid cells, and 12 distinct tumor cell states. LR scores were considered significant if p-value < 0.05, determined by 1000 cell label shuffles. LR scores were set to 0 for non-significant interactions and further averaged across samples to minimize sample bias. Next, LR scores between different tumor cell states were compared; pairs with a log2-fold change ≥ 0.5 and a signal in ≥ 10% of tumor cell states were selected. Finally, for each tumor cell state, a ranking score (score = P * FC) was calculated for each LR interaction based on the proportion of cells with the LR interaction (P) and the log2-fold change compared with other tumor cell states (FC), both rescaled between 0 and 1. The top three LR interactions for each tumor cell state were selected based on these ranking scores.

Molecular co-dependencies in tumor cell villages
We developed a bioinformatic method to determine the molecular co-dependencies between tumor cells and their neighboring non-tumor cells (including CAFs, TAMs, TECs, T cells, and B cells) within each tumor cell village. Specifically, for a given tumor-neighboring cell type pair in a tumor cell village, we identified the top 50 highly variable genes (HVGs) for each cell type. For each tumor cell-neighboring cell type pair, we randomly selected 1000 tumor cells, 1000 paired neighboring cells of a given cell type located based on a 40 µm distance. Pearson correlation coefficients were calculated between each pair of tumor-HVGs and non-tumor cell HVGs (each cell type separately). We repeated the random selections for 1000 times. The top 50 gene pairs with the highest correlations were selected for each tumor cell village. We further determined the correlations of these gene pairs within each tumor cell villages based on 40–80, 80–160, and 160–320 µm distance, and beyond each tumor cell village. The correlation scores were further linearly rescaled to a range between 0 and 1 within each village.

Random forest model to assess molecular co-dependency
We used a random forest model to evaluate the importance of the identified gene pairs in each tumor cell village. A random forest model was built using expression profiles of the identified gene pairs. For each tumor cell, the expression of the paired genes in a corresponding neighboring cell type was averaged based on a 40 µm distance. We used two-thirds of the tumor cells in each village as the training set, with default parameters from the R package randomForest (v4.7-1.2). Tumor cell villages were used as classifier labels. The remaining one-third of the tumor cells served as the test set to evaluate model performance. To examine the importance of the gene pairs, we spatially shuffled the neighbors of tumor cells to construct perturbed neighboring expression profiles and assess the model’s performance on these perturbed profiles. Additionally, we spatially shuffled the neighbors of SPINK1+ tumor cells to specifically evaluate the impact of the SPINK1-related gene pairs in Villages 3 and 6. The model-building and neighbor perturbation steps were repeated 100 times to mitigate potential random effects and ensure robustness of the results.

Spatial co-localization of SPINK1+ spots and spots positive for given genes
For each sample in the 10X Visium dataset, we identified SPINK1+ spots as those with the top 25% SPINK1 expression across all spots. Genes of interest (GOI, e.g., BGN, COL1A1, COL1A2, COL4A1, COL5A1, IGFBP7, LUM, MGP, TPT1, and VIM)-positive spots were similarly defined using the top 25% expression threshold. For each SPINK1+ spot, we used the “kNN” function from the R package dbscan (v1.2-0) to identify the nearest GOI-positive spots and calculate their distances. The distances for all SPINK1+ spots were averaged as the observed distance between SPINK1+ spots and GOI-positive spots for each sample. To assess significance, we generated 1000 random spot lists, each containing the same number of spots as GOI-positive spots. For each random list, “kNN” was used to calculate the average distance of SPINK1+ spots and the random spots (random distance). P-values were calculated as p = n/1000, where n is the number of random distances equal to or less than the observed distance for the given gene pair in a sample.

Gene set enrichment analysis (GSEA)
R package fgsea (v1.32.0) was used to perform GSEA analysis. Pathways with p-values < 0.05 were considered as significantly enriched pathways. For single sample GSEA (ssGSEA) analysis, R package GSVA (v2.3.1) was used. Pathways were obtained from the hallmark gene sets of GSEA database. Pathways with enrichment scores ≥ mean + standard deviation of all derived scores were retained.

Cell–cell communications based on single-cell transcriptome data
The function “nichenet_seuratobj_cluster_de” from the R package NicheNetR (v2.1.0) was used with default parameters to infer cell-cell communication in single-cell transcriptomics. For each of the identified ligands, the top-weighted receptor was selected from the NicheNet database. Ligand-receptor pairs were retained if the receptors (tumor cells as senders) or the ligands (CAFs as senders) were significantly differentially expressed between mCAF and other CAFs (adjusted p-value < 0.05). p-values for ligands were calculated using two-sided Wilcoxon test and corrected for multiple testing using the Benjamini–Hochberg false discovery rate.

Multiplex immunofluorescence
Snap-frozen CCA tissue was embedded in a tissue freezing medium (Leica) and cut into 5 µm sections. These sections were then stained using a standard immunofluorescence protocol. Primary antibodies against SPINK1 (1:50, Invitrogen, Cat # PA5-102353) and COL1A1 (1:80, Cell Signaling, Cat # 66948S) were used. The corresponding secondary antibodies conjugated with fluorophores (anti-rabbit Alexa 488, 1:500, Cat # 4412S; anti-mouse Alexa 555, Cell Signaling, Cat # 4409S) were then added. Counterstaining was performed with DAPI (1:1000, Carl Roth). The fluorescently stained tissues were viewed with a Leica DMI8 confocal microscope using a 63x magnification objective with the Tile Scan function, and the images were processed using LAS X Office software.

Statistics and reproducibility
No statistical method was used to predetermine sample size. No sample in this study was excluded from the analysis. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Liu et al. investigated the spatial organization of tumor cells and their microenvironment in primary liver cancer by integrating\nspatial single-cell imaging and single-cell RNA sequencing. They developed a deep learning-based strategy termed Spatial\nDynamics Network to spatially map tumor cell states and their surrounding architectures. With these new results, they\nintroduced the new concept of \"tumor cell villages\", which exhibited village-specific molecular co-dependencies between\ntumor cells and their microenvironment and were associated with patient outcomes. This concept highlights a centric role of\ntumor in tumorigenesis but also considers the critical contribution and diverse phenotype of immune and stromal cells.\nOverall, this paper is technically sound and potentially a great resource for liver cancer. The new tumor village concept is\ninteresting and offers a new angle to understand tumor heterogeneity. There are several concerns needed to be addressed\nbefore publication.\nMajor:\n1. The study includes both HCC and iCCA patients and notes specific enrichment patterns for tumor cell states and villages\nbetween these subtypes. Given the distinct biology and clinical management of HCC and iCCA, it would be beneficial to\nmore explicitly discuss whether certain villages or co-dependencies are predominantly found in one subtype versus the\nother and the potential clinical implications.\n2. The features of Villages 1 and 2, specifically a vascularized environment, the EMT state, and enrichment of iCAF, are\nwell-established indicators of poor tumor prognosis. Why do patients enriched for V1 and V2 show a better prognosis\ncompared to other patients (Figure 5E)?\n3. Figure S7B, derived from the primary CosMx SMI dataset, suggests that many of the identified tumor cell villages exhibit a\nlargely patient-specific distribution, often being prominent in only a single patient or a small subset of patients within the\ninitial cohort. However, when the derived village gene signatures are applied to validation cohorts (as indicated in analyses\nof 10X Visium data, Figure S7G/H), a different pattern seems to emerge where a majority of patients/samples appear to be\npredominantly composed of signatures related to Village 1, Village 2, and Village 8.\nThese two sets of observations do not seem to fully align. Could the authors please clarify this apparent inconsistency? Is\nthis discrepancy primarily attributable to the inherent differences and resolutions between the various technology platforms?\nAlternatively, could it suggest potential limitations in the stability, robustness, or generalizability of the identified village gene\nsignatures when extrapolated to more diverse patient datasets?\n4. A mIF staining of SPINK1+ tumor cells is needed to validate the existence and the surrounding TME architecture of these\nmalignant cells.\n5. The findings, like spatial preference of tumor cell state, were validated in published 10X Visium spatial transcriptome data\nfrom 39 HCC patients, which is comparable to the discovery cohort. However, the bulk RNAseq data used for validation\nlacks single cell resolution. A much larger single cell RNA-seq data may be more appropriate to validate the findings.\n6. How to explain the clinical relevance of the identified tumor villages since multiple villages are observed in one patient?\n7. The tumor niche or microenvironment is a well-established concept to describe the tumor cell and surrounding immune\ncells, etc. A special paragraph is needed to discuss the similarities and differences between the tumor niche and tumor\nvillage concepts to help the readers understand this proposed new concept.\nMinor:\n1. For the reader's convenience, within the same figure (for example, Figure 2C), the same cell type should be marked with\nthe same color. This will make it easier for the reader to match the cell's name with its position in the UMAP plot.\n2. Clustering of both tumor cells and non-tumor cells are not very clear to separate distinct cell types. There's a concern that\nbatch effects in the spatial single-cell data may have been over-corrected during processing, potentially affecting the\nclustering result.\n3. In the abstract section, it is the tumor village, not the Spatial Dynamics Network, which is obtained through the deep\nlearning-based strategy, as referred to the Methods sections.\n4. Typo: page 29, “recal” should be “recall”.", "reviewer_2": "In this manuscript, Liu et al have performed spatial transcriptomics and single cell RNA sequencing to interrogate the tumour\nheterogeneity, focussing on primary liver tumours. They documented that tumours are composed by distinct units constituted\nby different cell subtypes. In particular, different tumour cell state are preferentially associated with specific spatial dynamics\nnetworks. Collectively, results indicate a close relationship between tumour cell states and the cell subtypes composing their\nmicroenvironment. Outcomes propose a functional spatial context in which cell types composing each unit may mutually\ninfluence their support and behaviour. The authors introduced the concept of tumour cell “villages”, as distinct cell subtypes\nmay be linked to a reciprocal coexistence to support their entity, fitness, identify, and dependency. The functional relevance\nof these “villages” is exemplified by showing correlations between the presence of some villages and patient outcomes.\nMoreover, the authors exemplified a putative functional relationship by focussing on a gene, SPINK1, which expression in\ncancer cells of a given village appears to confer co-dependency with mCAFs.\nThe manuscript is very well structured, and the idea of tumour cell “villages” is very provocative and fascinating. Findings\ncontribute to functionally understand and interpret inter- and intra- tumour heterogeneity, its relevance in supporting distinct\ncancer cell subtypes, and a putative vulnerability in tumour heterogeneity, which could be exploited through strategies\ndesigned to perturb the strength of distinct “villages”.\nIn my opinion there is one major improvement needed to ensure clarity and avoid misinterpretations:\nSince the beginning, and already in the abstract, the authors speak about 50 samples. It is then going through the text,\nchecking details in supplementary information, that readers then realize that in fact these 50 samples correspond to 7\npatients. In the present version of the manuscript, this information is not evidently reported, as it appears only in Table S1.\nThe number of analysed patients must be clearly stated since the beginning. Moreover, the authors stated that they\nperformed studies on primary liver cancers. In fact, they used 3 iCCA and 4 HCC. iCCA and HCC have distinct molecular,\nhistological, vulnerable features. Moreover, male and female HCC have often very distinct characteristics. Consequently, the\nstudy is based on a small fraction of patients. Considering the complexity of the study and the deep analyses performed, I do\nnot find that 7 patients are not sufficient. The authors should clearly state these issues and clarify the samples on which the\nstudy has been performed, to avoid misunderstanding and misinterpretations from readers. I recommend as well to discuss\nthese limitations in the discussion sections, in addition to a pertinent section in which the authors have already discussed\nsome limitations of their study. I also want to stress that, despite the number and diversity of samples analysed could be a\nlimitation of the study, the main findings and concept of 'villages' are supported by their existence across different primary\nliver cancer, and notably in relation the complexity and heterogeneity of liver cancer. Moreover, the authors wrote in the\nresult section: “We noticed that tumor cell village compositions were much more similar across different sampling regions\nfrom the same patient compared to those from different patients, suggesting relatively stable tumor cell populations within\nindividual patients”. As I mentioned above, this may be again related to the comparison of patients with distinct primary liver\ncancer types (HCC and iCCA) and from male versus female patients. The authors should widen their interpretation, as the\none currently proposed would be only supported by analysing several patients belonging to the same primary liver cancer\ntype and sex of patients. This additional aspect could be as well elaborated in the discussion section.\nTo help readers interpreting differences identified in studies, I recommend to reports in all panels in figures and\nsupplementary figures, which samples belong to HCC and to iCCa patients.\nMinor comments:\n1) Genes should be in Italic.\n2) in S1, \"B\" should be positioned higher to make evident for readers that \"fraction of cells\" and \"average expression\" scales\nare referred to graphs below.", "reviewer_3": "In this manuscript, the authors collected a total of 50 samples containing tumor regions, boundary regions, and normal\nregions from tumor tissues of 3 iCCA patients and 4 HCC patients and studied intratumoral heterogeneity of primary liver\ncancer using CosMx and scRNA-seq technologies. First, the authors classified and annotated cancer cells and non-cancer\ncells. Second, the authors calculated the number of non-tumor cells within a 40 μm range around each tumor cell, used this\ninformation for cancer cell clustering, and defined the results as Spatial Dynamic Networks (SDNs). Combining cancer cell\nstates and SDNs, the authors constructed tumor villages using deep learning neural networks and found that Village 1 and\nVillage 2 contained mixed tumor cell states and vascularized microenvironments, which were associated with good\nprognosis. Finally, the authors explored the molecular co-dependencies maintaining tumor village stability. From my\nperspective, the intention of the study is meaningful, but some conclusions in the article are not solid, with specific issues as\nbelow.\nMajor concerns:\n1) In Figure 1A, the authors described multi-point tumor tissue sampling from the cancer region, boundary zone, and\nadjacent normal area for each patient. However, the paper provides minimal discussion on the differences among these\nthree regions and the spatial variations within the cancer region itself. An in-depth analysis of these differences might better\nreveal intratumoral heterogeneity and address actual biological and clinical problems. This would be more meaningful than\nmerely optimizing existing computational methods or showcasing fancier computational approaches.\n2) The tumor cell states obtained through NMF in Figure 2A contain multiple stress states. The first question is, given that\nStress/Metal and Stress are adjacent, why separate instead of combining them into a single Stress state? The second\nquestion is, the stress comprised of four clusters shows results in Figure S2E that are mainly concentrated in the 2CB\nsamples. Is this a sample-specific phenomenon? Including such information may not be conducive to explaining\nuniversal phenomena. The authors should filter out patient- or sample-specific states before proceeding with subsequent analysis.\n3) Can the numerous cell subtypes identified in Figure 2C be mapped to the cell subtypes currently recognized in the pan-\ncancer cell atlas?\n4) As shown in Figure 5C, the authors ultimately identified 8 tumor villages. In Figure S7B, the authors demonstrated the\ntumor village composition of each sample, showing that V1 is a village specific to patient 2H, while V2 is specific to patient\n3H. Looking at the SDNs in Figure 5C, T3, 4, and 5 are specific to V1, while T1 and 2 are specific to V2. Since I couldn't find\nthe composition of T across samples, I can only speculate whether V1 (T3, 4, 5) and V2 (T1, 2) respectively represent spatial\npatterns specific to patients 2H and 3H? If this is indeed the case, it seems unreasonable to combine these two patient-\nspecific tumor villages into one group for prognostic analysis.\n5) Furthermore, according to the description, both V1 and V2 belong to a vascularized environment. Theoretically, an\nenvironment associated with vascularization should be related to tumor metastasis and consequently linked to poor\nprognosis. However, the authors’ research results show that this is associated with good prognosis - how can this\nphenomenon be explained?\n6) The author introduces the concept of molecular co-dependencies at the end of the article, identifying multiple gene pairs\nthrough correlation analysis and suggesting these gene pairs are factors that maintain village stability. However, I notice that\nthese gene pairs are not ligand-receptor pairs. If the genes within these pairs cannot directly interact with each other, then\nthe claim that such pairing can stabilize the village seems rather far-fetched. This correlation is more likely due to spatial\nproximity, leading to related phenotypes in their respective cell types. Therefore, rather than attributing it to the action of gene\npairs, it would be more accurate to say that the coexistence of two cell types facilitates the formation of villages. I believe that\naccording to the concept proposed by the author, the gene pairs identified through analysis need to be explained through\nmolecular and cellular mechanisms, which currently seem weak and require more direct evidence.\nMinor concerns:\n1) Please label the sample size information from Table S1 in Figure 1A. While the article studies intratumor heterogeneity,\nthe schematic diagram creates confusion by showing only one sample each from tumor, boundary, and normal tissue.\n2) While I understand the author's intention to make the visualization more prominent, the scores in the six graphs in Figure\n1D have been adjusted too artificially and no longer appear to be continuous values.\n3) The authors present multiple clustering results in the paper, all simply named c1, c2, etc. Although Figure S4 indicates the\nsubtype classification of each cell type, Figure 2A does not clearly specify that these are tumor subtypes. It is suggested to\nadd abbreviations of the major categories before c to avoid confusion between different clustering results."}, "author_response": "Responses to Reviewers' Comments\nWe would like to express our sincere gratitude to the reviewers for their constructive and helpful\ncomments to improve our paper. Below are the point-by-point responses to the reviewers.\nReviewer #1:\nMajor comments:\n1. The study includes both HCC and iCCA patients and notes specific enrichment patterns for\ntumor cell states and villages between these subtypes. Given the distinct biology and clinical\nmanagement of HCC and iCCA, it would be beneficial to more explicitly discuss whether certain\nvillages or co-dependencies are predominantly found in one subtype versus the other and the\npotential clinical implications.\nResponse: Thank you for this helpful suggestion. We fully agree with the reviewer that HCC and\niCCA represent two clinically distinct subtypes of liver cancer, each with characteristic histological\nfeatures and clinical management strategies. At the same time, prior studies have demonstrated\nthat these subtypes may share certain molecular features (Chaisaingmongkol et al. Cancer Cell\n2017, PMID: 28648284). In line with the reviewer’s suggestion, we have expanded the Discussion\nto more explicitly elaborate on tumor cell villages in HCC and iCCA, the underlying molecular\nco-dependencies, and their potential clinical implications. In our analysis of single-cell spatial\ntranscriptomic data from the discovery cohort, we found that while most tumor cell villages were\ndetectable in both HCC and iCCA, certain villages may be preferentially enriched. Given the\nlimited number of patients in this cohort, we further examined village enrichment using a larger\nsingle-cell transcriptome dataset that included both HCC and iCCA patients. This analysis showed\nthat Villages 1 and 2 were preferentially enriched in HCC, whereas Villages 3–8 were observed in\nboth HCC and iCCA, as determined by hierarchical clustering of tumors using surrogate village\nfeatures. These results are consistent with prior evidence that HCC and iCCA have both shared\nand distinct features (Chaisaingmongkol et al. Cancer Cell 2017, PMID: 28648284; Ma et al. Nat\nCommun, 2022, PMID: 36476645). Survival analysis further revealed that patients with Villages\n1 and 2-enriched tumors had better overall survival outcomes than those enriched for Villages 3–\n8. Consistent results were observed when the analysis was restricted to HCC patients alone and\nfurther validated in three patient cohorts with bulk transcriptome data. We didn’t perform survival\nanalysis on iCCA alone in the single-cell cohort, as all iCCA samples in this dataset were enriched\nfor Villages 3–8. We also determined the molecular wiring underlying individual tumor cell\nvillages by identifying correlated gene pairs. Among them, critical correlations between SPINK1\nfrom tumor cells and genes from mCAFs were uncovered in Villages 3 and 6, which were linked\nto poor patient outcomes. This observation was further supported by cell–cell interaction analysis\nand validated using multiplex immunofluorescence staining. Moreover, in silico perturbation of\nSPINK1 and its paired genes resulted in significant decline of village identity prediction accuracy\ncompared with the original village configurations. These results highlight the important role of\n1\nthese gene pairs in maintaining the integrity of their corresponding tumor cell villages and pave\nthe way for developing new treatment strategies aiming at disrupting these co-dependencies to\ndestabilize the tumor. We have revised the manuscript accordingly. Please refer to Figures 5E-5H,\n6I, 6J, S11C-S11E, S12B, and the Discussion section.\n2. The features of Villages 1 and 2, specifically a vascularized environment, the EMT state, and\nenrichment of iCAF, are well-established indicators of poor tumor prognosis. Why do patients\nenriched for V1 and V2 show a better prognosis compared to other patients (Figure 5E)?\nResponse: We appreciate the reviewer for this insightful comment. To investigate this question, we\nperformed gene set enrichment analysis and further examined specific tumor and\nmicroenvironmental features linked to the identified tumor cell villages. We found that Villages 1\nand 2 were enriched in metabolic-related pathways, such as bile acid metabolism, adipogenesis,\nand fatty acid metabolism. In contrast, pathways such as MYC targets, E2F targets, G2/M\ncheckpoint, mitotic spindle, reactive oxygen species pathway, PI3K/AKT/mTOR signaling, and\nangiogenesis were observed across multiple Villages within 3–8, corresponding to programs that\npromote cell growth, survival, resistance to cell death, and neovascularization. Notably, while\nVillages 1 and 2 contained more endothelial cells, the endothelial cells in Villages 3–8 exhibited\nsignificantly higher expression of tumor-specific endothelial signatures (derived from Croix et al.\nScience 2000, PMID: 10947988), consistent with the angiogenesis pathway observed in several of\nthese villages. Previous single-cell studies have demonstrated that endothelial cells in cancer are\nheterogeneous, with some subtypes resembling normal vasculature and others displaying tumor-\npromoting properties (Li et al. Natl Sci Rev 2024, PMID: 39345334). Our results suggest that the\nendothelial cells in Villages 3–8 acquire more tumor-associated features compared with those in\nVillages 1 and 2, which may underlie the poorer prognosis observed in Villages 3–8. In our initial\nclassification of CAFs, we defined subtypes based on clustering analysis and the differentially\nexpressed genes in each cluster, following previously reported classifications of fibroblasts in the\nliterature (Cords et al. Nat Commun 2023, PMID: 37463917). Specifically, clusters enriched for\ninflammatory-related genes were annotated as iCAFs (CXCL12), rCAFs (CCL19, CCL21, CCL2),\nand ifnCAFs (IL32, CXCL9, CXCL10). After considering the reviewer’s comment, we revisited\nthe literature and found that the defined rCAFs and ifnCAFs in our study are generally considered\nas subpopulations within the broader iCAF category, which we had initially overlooked (Lavie et\nal. Nat Cancer 2022, PMID: 35883004). To improve clarity and avoid confusion, we have therefore\nrenamed these fibroblast subtypes according to their dominant marker genes: iCAFs as CXCL12+\nCAFs, ifnCAFs as IL32+ CAFs, and rCAFs as CCL19+ CAFs. The total abundance of the three\nCAF subtypes did not differ significantly between Villages 1–2 and Villages 3–8. With respect to\nEMT, although a small proportion of EMT-like cells were detected in Village 2, EMT pathways\nwere not significantly enriched, suggesting that these cells constitute only a minor fraction of the\ntumor and may not represent the dominant ecosystem state. These findings further underscore the\ncomplex interplay within the tumor ecosystem and highlight the value of the tumor cell village\n2\nconcept for understanding how coordinated tumor cell states and their surrounding environments\ncollectively shape overall tumor fitness and clinical outcomes. We thank the reviewer again for\nraising this question and we have revised the manuscript to include new analysis results. Please\nrefer to Figures 2C, S3, S11F-11H, as well as the corresponding revised text in the main\nmanuscript.\n3. Figure S7B, derived from the primary CosMx SMI dataset, suggests that many of the identified\ntumor cell villages exhibit a largely patient-specific distribution, often being prominent in only a\nsingle patient or a small subset of patients within the initial cohort. However, when the derived\nvillage gene signatures are applied to validation cohorts (as indicated in analyses of 10X Visium\ndata, Figure S7G/H), a different pattern seems to emerge where a majority of patients/samples\nappear to be predominantly composed of signatures related to Village 1, Village 2, and Village 8.\nThese two sets of observations do not seem to fully align. Could the authors please clarify this\napparent inconsistency? Is this discrepancy primarily attributable to the inherent differences and\nresolutions between the various technology platforms? Alternatively, could it suggest potential\nlimitations in the stability, robustness, or generalizability of the identified village gene signatures\nwhen extrapolated to more diverse patient datasets?\nResponse: We appreciate the reviewer for the comment. As noted by the reviewer, one contributing\nfactor is the resolution limitation of the 10x Visium platform, in which each spot measures 55 µm\nin diameter and contains a mixture of cell types. This inherent averaging effect reduces the ability\nto resolve fine-grained cell–cell relationships and accurately map tumor cell villages at true single-\ncell resolution, as achieved with CosMx SMI. Additionally, our initial 10x Visium analysis\nincluded only three HCC datasets available at the time, in which we observed enrichment of\nVillages 1, 2, and 7 (likely perceived by the reviewer as Village 8 due to color coding; we have\nmodified the colors to improve readability), consistent with our CosMx-based findings. To expand\nvalidation, we have now incorporated 10x Visium data from an iCCA cohort reported in a recent\npublication (Mo et al. Nature 2024, PMID: 39478210). Consistent with the CosMx results, we\nobserved that certain tumor cell villages, including Villages 3, 5, and 6, were enriched. These\nresults support the robustness of the identified tumor cell villages in our study. The new data are\npresented in Figure S11B.\n4. A mIF staining of SPINK1+ tumor cells is needed to validate the existence and the surrounding\nTME architecture of these malignant cells.\nResponse: Thank you for this helpful suggestion. Accordingly, we validated the colocalization of\nSPINK1+ tumor cells and mCAFs by performing co-staining of SPINK1 and COL1A1 (the most\nhighly expressed mCAF marker). Consistent with our spatial analysis, we observed strong\ncolocalization of SPINK1 and COL1A1 in tumors enriched for Village 6, where SPINK1+ tumor\ncells and mCAFs are closely associated. The data are presented in Figure S12B.\n3\n5. The findings, like spatial preference of tumor cell state, were validated in published 10X Visium\nspatial transcriptome data from 39 HCC patients, which is comparable to the discovery cohort.\nHowever, the bulk RNAseq data used for validation lacks single cell resolution. A much larger\nsingle cell RNA-seq data may be more appropriate to validate the findings.\nResponse: We appreciate the reviewer’s helpful suggestion and agree that bulk transcriptome data\noffers limited resolution. Accordingly, we performed additional validation of our findings using a\nsingle-cell RNA-seq dataset from 37 HCC or iCCA patients (Ma et al. J Hepatol 2021, PMID:\n34216724). Consistent with our initial results, patients enriched for Villages 3–8 had significantly\npoorer outcomes compared with those enriched for Villages 1 and 2. We further repeated the\nanalysis in HCC patients alone and observed a similar trend. We didn’t analyze iCCA separately\nbecause all iCCA tumors in this single-cell cohort were enriched for Villages 3–8. These new data\nare presented in Figures 5E, S11C, and S11D.\n6. How to explain the clinical relevance of the identified tumor villages since multiple villages are\nobserved in one patient?\nResponse: Thank you for this insightful question. As the reviewer noted, multiple villages can\ncoexist within a single tumor, reflecting the well-documented intratumor heterogeneity within\nindividual liver cancer patients. To assess the clinical relevance of this heterogeneity, we performed\nhierarchical clustering of tumors based on village-specific features. We found that despite the\ncoexistence of multiple villages within individual tumors, patients could be stratified into two\nmajor groups with significantly different outcomes: one enriched for Villages 1 and 2 features and\nthe other enriched for Villages 3–8 features. These findings suggest that while intratumor\nheterogeneity exists, certain dominant village-associated programs prevail at the patient level,\nallowing meaningful clinical stratification.\n7. The tumor niche or microenvironment is a well-established concept to describe the tumor cell\nand surrounding immune cells, etc. A special paragraph is needed to discuss the similarities and\ndifferences between the tumor niche and tumor village concepts to help the readers understand this\nproposed new concept.\nResponse: We thank the reviewer for this thoughtful suggestion. The concept of the tumor\nmicroenvironment or niche is well established. The tumor microenvironment describes the broader\necosystem of non-malignant components surrounding malignant cells, while tumor niche refers to\nmicroanatomical structures that maintain specific tumor cells, such as stemness, dormancy, or\nmetastatic potential. However, these concepts do not fully capture the coordinated spatial\norganization of diverse tumor cell states and their dynamic interactions with the\nmicroenvironment. To address this gap, we developed a graph attention neural network–based\n4\nmethod to identify spatially coordinated groups of tumor cell states, which we termed tumor cell\n“villages.” The village concept was inspired by human villages, which historically emerged to\nfulfill social, economic, and environmental needs, fostering community, adaptation, and resilience.\nAnalogously, tumor cell villages focus specifically on malignant cells as the central actors,\nrepresenting functional communities of diverse tumor cell states that are spatially organized and\nsupported by dynamic microenvironments. Heterogeneous tumor cell states within a tumor cell\nvillage may act collectively to enhance survival, adaptation, and progression, thereby reducing the\nvulnerability of individual cells. We have added a paragraph in the Discussion section accordingly\nto further elaborate on these concepts and to help readers better understand this proposed new\nconcept.\nMinor comments:\n1. For the reader's convenience, within the same figure (for example, Figure 2C), the same cell\ntype should be marked with the same color. This will make it easier for the reader to match the\ncell's name with its position in the UMAP plot.\nResponse: We thank the reviewer for this helpful comment. We have reviewed all relevant figures\nand standardized the color schemes so that the same cell type is represented by the same color\nwithin each figure to improve readability.\n2. Clustering of both tumor cells and non-tumor cells are not very clear to separate distinct cell\ntypes. There's a concern that batch effects in the spatial single-cell data may have been over-\ncorrected during processing, potentially affecting the clustering result.\nResponse: Thank you for this comment. To address the concern of potential over-correction, we\nperformed clustering of each non-malignant cell type for individual patients separately. We then\ncompared these patient-derived clusters with the clusters obtained from Harmony integration.\nOverall, we found that patient-derived clusters predominantly mapped to single Harmony clusters,\nindicating that Harmony preserved biologically meaningful groupings. For malignant cells,\nclustering was performed based on module scores of distinct cellular state-related gene programs\nderived from the original dataset, without additional Harmony integration. This approach, which\nhas been successfully applied in single-cell and spatial studies, allows us to identify shared\ntranscriptional programs across patients (Barkley et al. Nat Genet 2022, PMID: 35931863;\nGreenwald et al. Cell 2024, PMID: 38653236). We have added the new analysis results to Figure\nS4.\n3. In the abstract section, it is the tumor village, not the Spatial Dynamics Network, which is\nobtained through the deep learning-based strategy, as referred to the Methods sections.\nResponse: Thank you for this comment. We have modified it accordingly.\n5\n4. Typo: page 29, “recal” should be “recall”.\nResponse: Thank you for pointing out the typo. We have revised accordingly.\nReviewer #2:\n1. One major improvement needed to ensure clarity and avoid misinterpretations:\nSince the beginning, and already in the abstract, the authors speak about 50 samples. It is then\ngoing through the text, checking details in supplementary information, that readers then realize\nthat in fact these 50 samples correspond to 7 patients. In the present version of the manuscript, this\ninformation is not evidently reported, as it appears only in Table S1. The number of analysed\npatients must be clearly stated since the beginning. Moreover, the authors stated that they\nperformed studies on primary liver cancers. In fact, they used 3 iCCA and 4 HCC. iCCA and HCC\nhave distinct molecular, histological, vulnerable features. Moreover, male and female HCC have\noften very distinct characteristics. Consequently, the study is based on a small fraction of patients.\nConsidering the complexity of the study and the deep analyses performed, I do not find that 7\npatients are not sufficient. The authors should clearly state these issues and clarify the samples on\nwhich the study has been performed, to avoid misunderstanding and misinterpretations from\nreaders. I recommend as well to discuss these limitations in the discussion sections, in addition to\na pertinent section in which the authors have already discussed some limitations of their study. I\nalso want to stress that, despite the number and diversity of samples analysed could be a limitation\nof the study, the main findings and concept of 'villages' are supported by their existence across\ndifferent primary liver cancer, and notably in relation the complexity and heterogeneity of liver\ncancer. Moreover, the authors wrote in the result section: “We noticed that tumor cell village\ncompositions were much more similar across different sampling regions from the same patient\ncompared to those from different patients, suggesting relatively stable tumor cell populations\nwithin individual patients”. As I mentioned above, this may be again related to the comparison of\npatients with distinct primary liver cancer types (HCC and iCCA) and from male versus female\npatients. The authors should widen their interpretation, as the one currently proposed would be\nonly supported by analysing several patients belonging to the same primary liver cancer type and\nsex of patients. This additional aspect could be as well elaborated in the discussion section.\nTo help readers interpreting differences identified in studies, I recommend to reports in all panels\nin figures and supplementary figures, which samples belong to HCC and to iCCA patients.\nResponse: We sincerely appreciate the reviewer's thoughtful comments and valuable suggestions.\nWe also thank the reviewer for the positive assessment of our tumor cell village concept and the\ndata analyses supporting it. We fully agree that clarity regarding the number and types of patients\nanalyzed is critical to avoid misinterpretation. Accordingly, we have revised the figures and main\ntexts to clearly indicate the number of HCC and iCCA patients included in this study. In addition,\n6\nwe specify the clinical subtype of each sample in the main text and all figure panels, with sample\nIDs labeled with \"H\" for HCC and \"C\" for iCCA. We also acknowledge the limitation of our\nrelatively small discovery cohort, as well as the potential influence of factors such as patient sex,\netiology, and primary cancer type on tumor heterogeneity. Tumors from different patients evolve\nindependently, leading to unique gene expression profiles shaped by patient-specific\ncharacteristics. This heterogeneity between patients may reflect broader molecular, histological,\netiological, and demographic complexity, in addition to individual variations. These considerations\nfurther underscore the value of our tumor cell village concept. By spatially characterizing the\nlandscapes of individual tumors, the village concept enables a deeper understanding of how tumor\ncells establish unique landscapes within individual patients, while also facilitating comparisons of\nthese landscapes across patients. Moreover, defining tumor cell villages may uncover unique\nmolecular co-dependencies among cells that underlie the spatial organization of each village,\npotentially revealing new opportunities for therapeutic intervention. As the reviewer suggested,\nwe have expanded the Discussion section to more explicitly address the limitations of our study\nand we further emphasize the importance of understanding tumor spatial landscape through\ndefining tumor cell villages. Please refer to the revised Figure 1A and the revised text in the\nIntroduction, Results, and Discussion sections.\nMinor comments:\n1. Genes should be in Italic.\nResponse: Thank you for this comment. We have revised gene names in figures and the manuscript\naccordingly, using italic for gene names and regular font for proteins.\n2. in S1, \"B\" should be positioned higher to make evident for readers that \"fraction of cells\" and\n\"average expression\" scales are referred to graphs below.\nResponse: We have modified this figure accordingly. We thank the reviewer again for all the\nhelpful suggestions.\nReviewer #3:\nMajor comments:\n1. In Figure 1A, the authors described multi-point tumor tissue sampling from the cancer region,\nboundary zone, and adjacent normal area for each patient. However, the paper provides minimal\ndiscussion on the differences among these three regions and the spatial variations within the cancer\nregion itself. An in-depth analysis of these differences might better reveal intratumoral\nheterogeneity and address actual biological and clinical problems. This would be more meaningful\nthan merely optimizing existing computational methods or showcasing fancier computational\napproaches.\ndepth analysis to investigate spatial variations both among different tumor regions and within\nindividual regions. Specifically, we first conducted a detailed comparison of cellular compositions\nacross different tumor regions by examining the major cell types, malignant cell states, and non-\nmalignant cell subtypes. At the level of major cell types, we observed that T cells were markedly\ndecreased in the tumor border and core relative to the adjacent non-tumor region, whereas myeloid\ncells and fibroblasts progressively increased from the border to the core. Epithelial cells also\nshowed a gradual decrease in the border and core compared with the adjacent non-tumor region.\nDetailed comparison of tumor cell states revealed largely comparable proportions between the\nborder and core regions in our cohort. While most non-malignant cell subtypes showed no\nsignificant differences across these regions, several displayed region-specific enrichment.\nAt the tumor border, two immune subtypes (Kupffer cells and CD8+ effector T cells) and four stromal\nsubtypes (liver sinusoidal endothelial cells [LSECs], post-capillary venule endothelial cells,\nCCL19+ CAFs, and CXCL12+ CAFs) were significantly enriched compared to the tumor core.\nThis observation is consistent with the typical localization of Kupffer cells and LSECs in normal\nliver tissue. The additional enrichment of CD8+ effector T cells, inflammatory CAFs (CCL9+\nCAFs and CXCL12+ CAFs), and post-capillary venule endothelial cells highlights an inflamed\nmicroenvironment in this region. In contrast, tip cells and smooth muscle cells were elevated in\nthe tumor core, along with immune-regulatory TAMs (Reg-TAMs), suggesting a potential pro-\nangiogenic and immunosuppressive microenvironment.\nAdditionally, we investigated the spatial variation of cells within each tumor region by calculating\nMoran’s I to quantify spatial autocorrelation. For all major cell types, Moran’s I values were\npositive, indicating that each cell type exhibits some degree of spatial autocorrelation within\nindividual tumor regions. Notably, Moran’s I increased for epithelial cells, myeloid cells,\nfibroblasts, and endothelial cells in the tumor border and core compared with adjacent non-tumor\nregions, suggesting stronger spatial clustering of these populations within tumors. In contrast,\nMoran’s I values for T cells and B cells decreased in the tumor border and core relative to non-\ntumor regions, reflecting a more dispersed distribution within the tumor. These observations are\nconsistent with the clustered distribution of T cells and B cells around portal veins in normal liver,\nwhereas other cell types are more evenly dispersed in normal tissue. Malignant cell states displayed\nlow but positive Moran’s I values, with no significant differences between border and core,\nsuggesting only marginal spatial autocorrelation of each malignant cell state within individual\ntumor regions for the samples in our cohort. These new data have been shown in Figures S6 and\nS7 and we have also revised the texts in the main manuscript accordingly.\n2. The tumor cell states obtained through NMF in Figure 2A contain multiple stress states. The\nfirst question is, given that Stress/Metal and Stress are adjacent, why separate them instead of\ncombining them into a single Stress state? The second question is, the stress comprised of four\nclusters shows results in Figure S2E that are mainly concentrated in the 2CB samples. Is this a\nsample-specific phenomenon? Including such information may not be conducive to explaining\nuniversal phenomena. The authors should filter out patient- or sample-specific states before\nproceeding with subsequent analysis.\nResponse: Thank you for this comment. When determining malignant cell state, we calculated\nmodule scores for each tumor cluster and further annotated each cluster based on the dominant\nmodules with the highest scores. As shown by the log2 fold change of Stress and Metal module\nscores below, Cluster 10 had comparable enrichment for both modules, whereas Clusters 16, 21,\n24, and 7 showed much higher Stress scores than Metal scores. Therefore, we annotated Cluster\n10 as Stress/Metal and others as Stress.\nIn Figure S2E, the Stress state appears as a dominant malignant cell state only in sample 2CB.\nThis bar plot reflects the relative proportions of different cellular states within each sample. While it\nallows us to compare overall cellular compositions, it does not reflect the actual number of cells\nacross samples. To address this limitation, we have now provided a table reporting the absolute\nnumber of malignant cells in each state across individual samples. As shown in this table, Stress\nstate cells are present across different samples and their absolute numbers are relatively large,\nsuggesting that Stress state cells are not sample-specific. The data are presented in Table S4.\n3. Can the numerous cell subtypes identified in Figure 2C be mapped to the cell subtypes currently\nrecognized in the pan-cancer cell atlas?\nResponse: Thank you for raising this important question. We annotated cell subtypes within each\nmajor non-malignant cell type based on differentially expressed genes and further used publicly\navailable cell atlases as references. Accordingly, the identified cell subtypes are well mapped to\npublicly available pan-cancer cell atlases. Specifically, for endothelial cell (EC) subtypes, tip cells\n(CXCR4/collagen genes), arterial ECs (EFNB2), and lymphatic ECs (PROX1) have been identified\n9\nin across multiple cancer types (Li et al. Natl Sci Rev 2024, PMID: 39345334; Goveia et al. Cancer\nCell 2020, PMID: 31935371; Vanlandewijck et al. Nat 2018, PMID: 29443965; Zhao et al. Nat\nCancer 2024, PMID: 38172338). Postcapillary veins (VCAM1/CLU), activated capillaries\n(CAV1/IGFBP7), CD14+ capillaries (CD14), and scavenging capillaries (CD68) have been\ndescribed in lung cancer (Goveia et al. Cancer Cell 2020, PMID: 31935371). LSEC (LYVE1) are\nunique to the liver (MacParland et al. Nat Commun 2018, PMID: 30348985). We also identified a\ncluster of endothelial cells showing high expression of RGCC but lacking well-established\nmarkers, and therefore we annotated the cluster as unknown. We identified eight subtypes within\nthe fibroblast population. Among them, mCAFs (collagen genes), apCAFs (CD74 and antigen-\npresentation), and three inflammatory subtypes (CCL9+ CAFs, IL32+ CAFs, and CXCL12+\nCAFs) expressing high level of inflammatory and chemo-attractive mediators have been described\nin various cancer types (Liu et al. Cancer Cell 2025, PMID: 40154487; Cords et al. Nat Commun\n2023, PMID: 37463917; Zhang et al. J Hepatol 2020, PMID: 32505533; Lavie et al. Nat Cancer\n2022, PMID: 35883004). eCAFs (KRT17) has been reported in iCCA (Zhang et al. J Hepatol 2020,\nPMID: 32505533), while tCAFs (VEGFA) have been identified in breast and lung cancers (Cords\net al. Nat Commun 2023, PMID: 37463917; Cords et al. Cancer Cell 2024, PMID: 38242124). In\naddition, smooth muscle cells (RGS5/ACTA2) are commonly observed in various normal and\ncancer tissues. Among the ten myeloid cell subtypes, nine subtypes including Prolif-TAM\n(HMGN2/TUBB), Reg-TAM (GPNMB), Angio-TAM (TIMP1), INT monocyte (SERPINA1), LA-\nTAM (APOE), TIM (CD14), and Inflam-TAM (CCL20), neutrophils (CSF3R) and mast cells\n(CPA3) have been reported in various cancer studies (Cheng et al. Cell 2021, PMID: 33545035,\nXue et al. Nat 2022, PMID: 36352227, Ma et al. Trends Immunol 2022, PMID: 35690521).\nKupffer cells (CD5L) are unique to the liver (Aizarani et al. Nat 2019, PMID: 31292543). For the\nT cell subtypes, T (GZMK/GZMB/GNLY), T (ANXA1), T (FOXP3), T (IL7R/GZMK),\neff mem reg EMRA\nProlif-T (MKI67) and IFN-I responsive T (IFIT1) were annotated based on gene markers obtained\nfrom a pan-caner study of T cells across 21 cancer types as well as T cell studies of individual\ncancer types (Zheng et al. Science 2021, PMID: 34914499; Zheng et al. Cell 2017, PMID:\n28622514; Ma et al. J Hepatol 2021, PMID: 34216724). LGALS1+, EOMES+, CXCR6+, and JUN+\nwere annotated based on their highly expressed genes. We also observed three clusters of T cells\nlacking classical markers and we therefore annotated them as unknown. The five B cell subtypes\nwere annotated based on the markers from pan-cancer studies of B cells (Yang et al. Cell 2024,\nPMID: 39047727; Ma et al. Science 2024 PMID: 38696569).\n4. As shown in Figure 5C, the authors ultimately identified 8 tumor villages. In Figure S7B, the\nauthors demonstrated the tumor village composition of each sample, showing that V1 is a village\nspecific to patient 2H, while V2 is specific to patient 3H. Looking at the SDNs in Figure 5C, T3,\n4, and 5 are specific to V1, while T1 and 2 are specific to V2. Since I couldn't find the composition\nof T across samples, I can only speculate whether V1 (T3, 4, 5) and V2 (T1, 2) respectively\nrepresent spatial patterns specific to patients 2H and 3H? If this is indeed the case, it seems\n10\nunreasonable to combine these two patient-specific tumor villages into one group for prognostic\nanalysis.\nResponse: We thank the reviewer for this comment. We agree with the reviewer that patient 2H is\ndominated by Village 1 and patient 3H by Village 2. However, we also observed that Villages 1\nand 2 were present at high proportions in other patients (Village 1 in patient 2C and Village 2 in\npatient 1H). We recognized that the relatively small number of patients in our discovery cohort is\na limitation, which may give the appearance that certain villages are restricted to only a few\npatients. To address this concern, we extended our analysis by applying Villages 1 and 2 signatures\nto independent datasets, including 10x Visium spatial transcriptome data from 46 samples across\nfour liver cancer studies (three HCC studies and one iCCA study), scRNA-seq data from 37 liver\ncancer patients, and bulk transcriptome data from 674 liver cancer patients. Across these larger\ndatasets, Villages 1 and 2 signatures were frequently observed across different patients,\ndemonstrating that these tumor cell villages are not specific to a few patients. Additionally, their\nconsistent association with clinical outcomes across different patient cohorts further supports that\nVillages 1 and 2 are biologically meaningful. The data are presented in Figures 5E-5H and S11A-\nS11E. We have revised the manuscript accordingly and also expanded the Discussion section to\nreflect the limitations of our study.\n5. Furthermore, according to the description, both V1 and V2 belong to a vascularized\nenvironment. Theoretically, an environment associated with vascularization should be related to\ntumor metastasis and consequently linked to poor prognosis. However, the authors’ research results\nshow that this is associated with good prognosis - how can this phenomenon be explained?\nResponse: We appreciate the reviewer's comment and agree with the reviewer that tumor\nangiogenesis is usually linked to poor patient outcomes. To investigate this question, we further\nexamined tumor vasculature in the identified tumor cell villages. We found that although Villages\n1 and 2 showed overall enrichment in endothelial cells, the expression of tumor-specific\nendothelial signature genes (derived from Croix et al. Science 2000, PMID: 10947988) was\nsignificantly higher in Villages 3–8 compared with Villages 1 and 2. Previous single-cell studies\nhave demonstrated that endothelial cells in cancer are heterogeneous, with some subtypes\nresembling normal vasculature and others displaying tumor-promoting properties (Li et al. Natl\nSci Rev 2024, PMID: 39345334). Our results suggest that the endothelial cells in Villages 3–8 may\nacquire more tumor-associated features, which may underlie the poorer prognosis observed in\nthese patients. The new data are presented in Figure S11G.\n6. The author introduces the concept of molecular co-dependencies at the end of the article,\nidentifying multiple gene pairs through correlation analysis and suggesting these gene pairs are\nfactors that maintain village stability. However, I notice that these gene pairs are not ligand-\nreceptor pairs. If the genes within these pairs cannot directly interact with each other, then the\n11\nclaim that such pairing can stabilize the village seems rather far-fetched. This correlation is more\nlikely due to spatial proximity, leading to related phenotypes in their respective cell types.\nTherefore, rather than attributing it to the action of gene pairs, it would be more accurate to say\nthat the coexistence of two cell types facilitates the formation of villages. I believe that according\nto the concept proposed by the author, the gene pairs identified through analysis need to be\nexplained through molecular and cellular mechanisms, which currently seem weak and require\nmore direct evidence.\nResponse: We greatly appreciate the reviewer's insightful comment. Our initial motivation was to\nunderstand molecular co-dependencies within individual tumor cell villages that extend beyond\nligand–receptor interactions. While ligand–receptor signaling captures an important subset of\nintercellular communication, cells also interact through additional mechanisms such as\nextracellular vesicles, metabolites, adhesion molecules, gap junctions, and nanotubes. By\nexamining molecular co-dependencies, our goal was to capture broader molecular relationships\nthat may underlie village stability. We fully agree with the reviewer that the gene pairs we\nidentified were based on spatial correlation and that elucidating the underlying molecular and\ncellular mechanisms is critical for understanding the biology of tumor cell villages. In the original\nmanuscript, we highlighted SPINK1 in tumor cells and its associated genes in mCAFs as a proof\nof concept. We further demonstrated that SPINK1+ tumor cells engaged more extensively with\nmCAFs than with other CAF subtypes through ligand-receptor interactions. We provided only the\ntop 20 ligand–receptor interactions between SPINK1+ tumor cells and CAFs in the manuscript,\nleaving many additional interactions unexplored. In response to the reviewer's comment, we now\nprovide all the identified ligand–receptor interactions. Notably, SPINK1 itself was identified as\none of the ligands in these interactions with EGFR serving as its receptor. Consistently, SPINK1\ndownstream genes were expressed at significantly higher levels in mCAFs compared with other\nCAFs. Similarly, the EGFR downstream genes were significantly upregulated in mCAFs. We also\nperformed immunofluorescence staining for SPINK1 and COL1A1 (a representative marker for\nmCAF) to confirm the colocalization of the two proteins. Previous studies have demonstrated that\nSPINK1 can stimulate the proliferation of both tumor cells and fibroblasts through binding to\nEGFR and activating its downstream signaling using experimental approaches (Ateeq et al. Sci\nTransl Med 2011, PMID: 21368222; Ozaki et al, Mol Cancer Res 2009, PMID: 19737965). It is\ntherefore likely that SPINK1 may also promote the proliferation of mCAFs via EGFR signaling,\nproviding molecular mechanisms for the SPINK1-related gene pairs identified in our study. The\nnew results have been added to Figure S12 and Tables S7 and S8. As we only demonstrated the\nmolecular mechanisms of SPINK1 as a proof of concept in this study, we have expanded the\nDiscussion to acknowledge the limitation of our study and to emphasize the need for further\nexperimental approaches to elucidate the molecular and cellular mechanisms of the identified gene\npairs. We thank the reviewer again for this valuable comment, which is very helpful to improve\nour manuscript.\n1. Please label the sample size information from Table S1 in Figure 1A. While the article studies\nintratumor heterogeneity, the schematic diagram creates confusion by showing only one sample\neach from tumor, boundary, and normal tissue.\nResponse: Thank you for pointing this out. We have revised Figure 1A accordingly.\n2. While I understand the author's intention to make the visualization more prominent, the scores\nin the six graphs in Figure 1D have been adjusted too artificially and no longer appear to be\ncontinuous values.\nResponse: We appreciate the comment and have modified the color in Figure 1D to continuous\nvalues accordingly.\n3. The authors present multiple clustering results in the paper, all simply named c1, c2, etc.\nAlthough Figure S4 indicates the subtype classification of each cell type, Figure 2A does not\nclearly specify that these are tumor subtypes. It is suggested to add abbreviations of the major\ncategories before c to avoid confusion between different clustering results.\nResponse: Thank you for this comment. We have modified the cluster naming scheme in Figure\n2A, assigning the prefix “tc” to tumor clusters for clarity. We thank the reviewer again for all the\nconstructive and helpful comments to improve our paper.\n13"}, {"round": 2, "reviewer_comments": {"reviewer_1": "All my concerns have been addressed. I have no further comments.", "reviewer_2": "The authors have adequately addressed the points raised. Congratulation for this inspiring work.", "reviewer_3": "The authors have basically answered all the questions. However, I still have a minor question about the response to the first comment. Regarding the epithelial cells in the Supplementary Figure 6A, according to the description in line 579 of the manuscript, these cells include hepatocytes, cholangiocytes, and tumor cells. Among them, normal epithelial cells should be enriched in the adjacent non-tumor area, while tumor cells should be enriched in the tumor area. Therefore, observing the regional enrichment pattern of total epithelial cells, which include both normal and cancerous epithelial cells, seems unreasonable."}, "author_response": "Responses to Reviewers' Comments\nReviewer #1:\nAll my concerns have been addressed. I have no further comments.\nReviewer #2:\nThe authors have adequately addressed the points raised. Congratulation for this inspiring work.\nResponse: We thank the reviewers again for all the constructive and helpful comments to improve\nour paper.\nReviewer #3:\nThe authors have basically answered all the questions. However, I still have a minor question about\nthe response to the first comment. Regarding the epithelial cells in the Supplementary Figure 6A,\naccording to the description in line 579 of the manuscript, these cells include hepatocytes,\ncholangiocytes, and tumor cells. Among them, normal epithelial cells should be enriched in the\nadjacent non-tumor area, while tumor cells should be enriched in the tumor area. Therefore,\nobserving the regional enrichment pattern of total epithelial cells, which include both normal and\ncancerous epithelial cells, seems unreasonable.\nResponse: Thank you for this comment. We initially presented the total epithelial abundance in\nSupplementary Figure 6A to provide an overview of epithelial cell changes across different tumor\nregions. In response to the reviewer’s comment, we have now added the detailed abundances of\nmalignant and non-malignant (hepatocytes and cholangiocytes) epithelial cells. As expected, non-\nmalignant epithelial cells were predominantly located in the adjacent non-tumor region, whereas\nmalignant cells were mainly enriched in the tumor core. Please refer to the revised Supplementary\nFigure 6A. We thank the reviewer again for this helpful suggestion.\n1"}, {"round": 3, "reviewer_comments": {"reviewer_3": "I have no further comments."}, "author_response": "Responses to Reviewers' Comments\nReviewer #3:\nI have no further comments.\nResponse: We thank the reviewer again for all the helpful comments to improve our paper.\n1"}]

Decision letter