13.3
Impact Factor

PDF

Theranostics 2026; 16(7):3870-3891. doi:10.7150/thno.125381 This issue Cite

Research Paper

Pathogenomic analysis reveals clinically relevant epithelial-mesenchymal plasticity in esophageal squamous cell carcinoma

Ruzhen Chen^1,2,5#, Chenyi Xie^2,5#, Ziyu Ning⁶, Meng Yang³, Zezhuo Su⁴, Jiahui Chen², Kunheng Du², Yihuai Hu⁷, Chu Han^2,5, Shaojun Zhang⁸, Qingling Zhang⁷, Meng Liu³ , Zaiyi Liu^1,2,5

1. School of Medicine, South China University of Technology, Guangzhou 510006, China.
2. Department of Radiology, Guangdong Provincial People's Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou 510080, China.
3. Affiliated Tumor Hospital of Xinjiang Medical University, Urumqi 830011, Xinjiang, China.
4. Department of Orthopaedics and Traumatology, School of Clinical Medicine, Li Ka Shing Faculty of Medicine, The University of Hong Kong, Pok Fu Lam, Hong Kong SAR, China.
5. Guangdong Provincial Key Laboratory of Artificial Intelligence in Medical Image Analysis and Application, Guangzhou 510080, China.
6. Department of Pathology, Guangdong Provincial People's Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou 510080, China.
7. Department of Thoracic Surgery, Guangdong Provincial People's Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou 510080, China.
8. Guangdong Academy of Medical Sciences and Medical Research Institute, Guangdong Provincial People's Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou 510080, China.
^# Contributed equally to this work.

✉ Corresponding authors: Meng Liu, MD, PhD (liuxiaomeng198154com). Zaiyi Liu, MD, PhD (liuzaiyiorg.cn).

Received 2025-9-17; Accepted 2026-1-5; Published 2026-1-14

Abstract

Rationale: Esophageal squamous cell carcinoma (ESCC) is a highly aggressive malignancy. The metastasis and poor prognosis of ESCC are closely associated with tumor microenvironment (TME) heterogeneity, which is driven by epithelial-mesenchymal transition (EMT). Clinically, how to diagnose and target EMT progression remains a key challenge for ESCC.

Methods: Integration of pathological images and bulk RNA sequencing profiles identified a high-risk subtype exhibiting EMT enrichment and immunosuppression. Single-cell and spatial transcriptomics revealed EMT macrostates and their spatial distribution. The role of CACNA1C in programming malignant phenotype was tested in vitro. A pathological image-based deep learning model successfully predicted the spatial expression distribution of CACNA1C, indicating possible clinical utility.

Results: EMT progression comprised three macrostates: the early state (high epithelial and metastatic potential), the stable state (hybrid E/M phenotype and high stemness), and the late state (high mesenchymal and invasive propensity). ITGA3 and ITGB4 antagonistically regulate malignant phenotype in the early state. Notably, suppression of CACNA1C induced transdifferentiation from stable/late-state cells to normal epithelium-like cells.

Conclusions: This study provides novel insights into the EMT mechanism in ESCC, proposes an intervention strategy, and emphasizes the promising clinical application of pathological images in EMT assessment.

Keywords: esophageal squamous carcinoma, epithelial mesenchymal transition, histopathology, bulk sequencing, single-cell sequencing, spatial transcriptomes, deep learning

Introduction

Esophageal cancer is the seventh leading cause of cancer-related death worldwide [1]. Esophageal squamous cell carcinoma (ESCC) accounts for approximately 90% of esophageal cancer cases, and more than 70% of ESCC cases are locally advanced at the time of diagnosis [2]. ESCC is characterized by aggressive submucosal spread and lymphatic metastasis, and significantly poorer survival outcomes are observed in patients with metastatic disease [3]. While surgical resection is the fundamental treatment for ESCC [4], emerging therapeutic targets and molecular classifications derived from multiomics technologies have been increasingly proposed to guide personalized therapy [5]. Understanding the ESCC tumor microenvironment (TME) is essential for overcoming treatment resistance, refining subtype-specific treatment strategies, and improving patient outcomes.

The degree of epithelial-to-mesenchymal transition (EMT) has been reported to be a predictor of tumor invasion and lymph node metastasis [6, 7]. EMT is a cellular process in which polarized epithelial cells undergo extensive molecular reprogramming and phenotypic transition [8]. EMT occurs not only during embryonic development, tissue regeneration, and wound healing but also in tumors [9], which leads to a reduction in adhesion between cells and an increase in the ability to migrate, thus initiating the transformation from a benign tumor to an invasive tumor. Recent studies have shown that EMT encompasses a continuum of phenotypes, with multiple mixed states driving the transition from a complete epithelial phenotype to a complete mesenchymal phenotype [10, 11]. Elucidation of EMT differentiation trajectories is crucial for understanding ESCC invasiveness and metastasis and had significant potential to inform targeted therapies and improve patient survival; however, the specific EMT intermediate states and their functions remain incompletely defined.

Tumor morphology and architecture reflect biologically relevant subtypes, and integrating spatial pathology with genomic data enables deeper insights into the TME and biomarker development [12]. For example, Wang et al. and Yoo et al. performed unsupervised clustering on quantitative features characterizing the cell nuclei and the spatial relationships between cancer cells and lymphocytes to develop histological subtypes in liver cancer and colorectal cancer, respectively [13, 14]. Zhao et al. further incorporated stromal cell features and identified clinically relevant subtypes of breast cancer [15]. Recent deep learning advances have enabled cost-effective, reliable biological stratification directly from histopathological images in both clinical workflows and retrospective cohorts, revealing novel insights into tumor heterogeneity [16, 17]. The integration of pathogenomics and deep learning for multimodal biomarker discovery enables the characterization of tumor heterogeneity and clarification of previously unexplored genomic mechanisms in ESCC histological subtypes, advancing both oncology research and clinical decision-making.

In this study, we identified two pathological subtypes of ESCC with prognostic differences. We performed integrated analysis of genomics data for paired samples to comprehensively characterize the molecular dynamics and TME landscapes underlying individual subtypes, identified EMT as a critical determinant of prognosis, and mined therapeutic targets and functional ligand-receptor pairs for the high-risk subtype. In addition, we constructed ESCC-specific EMT trajectories and identified three macrostates with distinct functions. We found that ITGA3 and ITGB4 are functionally mutually exclusive and that CACNA1C is a promising trans-differentiation target for EMT. Finally, we developed a deep learning model to predict the expression and spatial distribution of CACNA1C. In conclusion, our study reveals the widespread impact of EMT on ESCC invasion, metastasis, prognosis, and treatment.

Methods

Multi-center data collection

This study has been approved by the institutional review boards at each hospital and informed consent was waived due to retrospective nature. This retrospective study involved a discovery cohort (Affiliated Tumor Hospital of Xinjiang Medical University, XJ, n = 125), as well as a validation cohort (Guangdong Provincial People's Hospital, GDPH, n = 78). Eligible participants met the following criteria: (i) a confirmed pathological diagnosis of ESCC; (ii) receipt of primary curative-intent surgical resection; and (iii) availability of comprehensive histological data. Patients were excluded if they had: (i) undergone any preoperative anti-cancer treatment; (ii) a prior history of malignancy or coexisting cancers; or (iii) missing or incomplete clinical, imaging, or pathological records. The overall survival (OS) as the primary endpoint, defined as the time interval from surgical resection to death from any cause. The minimum follow-up period was 36 months after surgery. We collected clinicopathologic characteristics from medical records and assessed followed-up information by pathology/imaging reports or telephone follow-up. Baseline characteristics were collected, including age, gender, tumor location, stage, T stage, N stage, and adjuvant therapy record. Statistical comparisons across cohorts were calculated by Chi-squared test (Table S1).

Public data collection

The published clinical specimens here were collected from the TCGA Research Network (https://portal.gdc.cancer.gov/) and the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). From the GDC Data Portal, we identified esophageal squamous cell carcinoma (ESCC) samples using the following selection criteria: PrimarySite = "Esophagus", DiseaseType = "Squamous cell neoplasms", SampleType = "Primary Tumor", and TreatmentOrTherapy = "No". We further refined our selection to include only cases with available H&E-stained diagnostic slides from surgical specimens. We employed the GDC Data Transfer Tool for download of clinical information and histopathological slides, while utilizing the TCGAbiolinks R package [18] to retrieve single nucleotide variation (SNV) data, copy number variation (CNV) data, and transcriptomic profiles. Additionally, we acquired single-cell RNA sequencing data (GSE160269) of ESCC clinical specimens from the GEO database.

Cell line data were obtained from multiple public resources: transcriptomic profiles from the Cancer Cell Line Encyclopedia (CCLE) [19] (https://sites.broadinstitute.org/ccle/), drug response data from the Genomics of Drug Sensitivity in Cancer (GDSC) [20] (https://www.cancerrxgene.org/), metastatic potential data from MetMap500 [21] (https://depmap.org/metmap/), genetic perturbation data from the Library of Integrated Network-based Cellular Signatures (LINCs) [22] (https://lincsproject.org/), and genetic dependency data from DepMap [23] (https://depmap.org/portal/data_page/).

H&E-stained section scanning

The formalin-fixed and paraffin-embedded tissue sections were collected by surgical resection and stained with H&E. A representative section of each patient was selected by experienced pathologists (Z.N. and Q.Z.) blinding of patient outcomes. The H&E-stained sections were scanned and digitalized using whole-slide scanners at 400×magnification with a resolution of 0.25 μm per pixel. All WSIs were manually checked for quality control, excluding images with blurry areas, and light- or over-stained areas.

Spatial transcriptomic sequencing

For 10x Visium v2 experiment, FFPE tissues were sectioned (5 μm), and following RNA quality assessment, sections underwent deparaffinization, H&E staining, and imaging. Probe hybridization employed the RTL technique, and ligation products were captured on Visium slides via poly(dT) sequences. Subsequent library construction included mRNA digestion, probe extension for barcode/UMI incorporation, alkaline elution, pre-amplification, indexing, and cleanup. Final libraries were sequenced on an Illumina platform.

For 10x Visium HD experiment, fresh tissue was fixed, paraffin-embedded, and sectioned. After deparaffinization, H&E staining, imaging, and decrosslinking per the manufacturer's protocol, probe hybridization and library construction were performed using the Visium HD kit. Libraries were sequenced in PE-150 mode on an Illumina platform.

For Stereo-seq experiment, fresh tissue was OCT-embedded, frozen, and cryosectioned (10 μm). Sections were adhered to a Stereo-seq poly-T chip, then processed for methanol fixation, H&E staining, imaging, permeabilization, reverse transcription, and cDNA purification according to kit instructions. Libraries were prepared and sequenced on a DNBSEQ-T7 instrument.

For 10x Xenium experiment, FFPE sections (5 μm) were mounted on Xenium slides. After histopathological assessment, spatial profiling was conducted using the Xenium platform with a 5K gene panel. The workflow included deparaffinization, probe hybridization, ligation, rolling circle amplification. Image fields of view were computationally stitched using the DAPI channel to reconstruct whole-tissue spatial maps. Cell segmentation was performed using a multimodal approach integrating nuclear and cytoplasmic signals, and detected transcripts were assigned to individual cells based on their spatial localization. Following completion of imaging, slides were subjected to post-run H&E staining and whole-slide scanning. H&E images were aligned with Xenium morphology images using Xenium Explorer software to enable visualization of spatial gene expression within the histological context.

Cell culture and transfection

Human ESCC cell line TE-8 (Ethephon, YCL-0533) was selected because of intrinsic CACNA1C expression allowing for knock-down of function assay. TE-8 was cultured in RPMI1640 media supplemented with 10% FBS and 1% P/S under 5% CO₂ conditions. Short tandem repeat profiling of the cell lines was performed to ensure integrity of cell line.

For knock-down assay, TE-8 was transfected with siRNA targeting CACNA1C or control siRNA (Hycyte, HX-H-S3-16863). At 48 hours post-transfection, RT-qPCR was performed to confirm that CACNA1C expression was decreased.

Cellular invasion assay

For cell invasion assay, cell lines were plated in a 24-well plate at 9.0 × 10⁴ cells per well. Invasion was evaluated using the Matrigel 24-well Cell Invasion Chamber (8 μm, Corning, 354480). After incubation at 37 °C with 5% CO₂ for 24 h, non-invasion cells on the upper membrane surface were removed by gentle scraping with a cotton swab. Membranes were then fixed in 4% formaldehyde (15 min, room temperature) and stained with 0.1% crystal violet (20 min). Invasive cells on the lower membrane surface were visualized under bright-field microscopy. Five different 200× filed images were taken for each well, and the number of invasive cells was counted by ImageJ (National Institutes of Health).

Cellular proliferation and colony formation assay

For cell proliferation assay, cell lines were plated in a 96-well plate at 1.0 × 10³ cells per well. Proliferation was evaluated using the CellTiter-Glo Luminescent Cell Viability Assay (Promega, G7570). After incubation at 37 °C with 5% CO₂ for 48 h, media was replaced to 200 μL of media containing 50% CellTiter-Glo reagent. Contents were mixed on an orbital shaker for 2 min, and incubated at room temperature for 10 min. The intensity of luminescent signal was detected by Luminometer.

For cell colony formation assay, cell lines were plated in triplicate in a 6-well culture plates at 1.0 × 10³ cells per well. Seven days after plating, cells were fixed and stained with crystal violet. The colony accounts were counted by ImageJ (National Institutes of Health).

Pathological feature extraction and analysis

The analysis pipeline began by segmenting whole-slide images (WSIs) into 4096 * 4096-pixel patches representing TME. Each patch was processed at 40× magnification using Hovernet [24], a pre-trained cell segmentation model on Panuke dataset [25], to identify tumor cells, immune cells, and stromal cells. Spatial positions of these cell types were extracted and used to construct topological graphs via minimum spanning trees. From these graphs, we derived six distinct cellular interaction subnetworks: tumor-tumor (T-T), tumor-lymphocyte (T-I), tumor-stroma (T-S), lymphocyte-lymphocyte (I-I), lymphocyte-stroma (I-S), and stroma-stroma (S-S). Patch-level features were aggregated into WSI-level representations using four statistical measures (mean, variance, skewness, and kurtosis), generating a final pathological feature matrix of 1,008 dimensions per WSI.

To evaluate cross-cohort feature consistency, we first applied the ComBat function from the sva R package [26] for batch effect correction, followed by PCA dimensionality reduction for visualization using the first two principal components. Intra-cohort heterogeneity was assessed by selecting the top 10 principal components contributing most to cohort variation (determined by the elbow method) and performing K-means clustering to stratify cases into two pathological subtypes. The clinical relevance of these phenotypes was evaluated through Kaplan-Meier survival analysis with log-rank tests, supplemented by Cox proportional hazards modeling to quantify survival impact, using p < 0.05 as the significance threshold.

Bulk transcriptomic feature extraction and functional profiling

Using the low-risk subtype as reference, we identified differentially expressed genes (DEGs) between pathological subtypes using the DESeq2 R package [27], obtaining log2 fold-change (logFC), p-values for each gene. BH method was applied for p-value significance correction. We performed gene set over representation analysis (ORA) on GO (Gene Ontology) terms from MSigDB, using significantly upregulated genes as the input gene list, implemented using the ClusterProfiler R package [28]. We conducted variation and enrichment analyses on gene sets from MSigDB, including Hallmark, KEGG, Reactome, CGP, and PID collections, implemented using the GSVA [29] and clusterProfiler R packages, respectively.

Pathology-pathway association

We employed regularized generalized linear models to investigate associations between pathological features and molecular pathways. Specifically, the model incorporated 1,008 WSI-level pathological features as predictors and pathway GSVA scores as response variables, implemented using the glmnet R package [30]. Model performance was evaluated via prediction accuracy in five-fold cross-validation, with lambda (regularization parameter) selection optimized through nested ten-fold cross-validation within each fold. Finally, we assessed model accuracy by computing Spearman's rank correlation coefficients between predicted pathway scores and actual GSVA scores for each pathway.

Disease-drug-target association

We computed drug-disease associations using the summarized Reverse Gene Expression Score (sRGES) established by Chen et al. [31] To validate the association reliability, we correlated sRGES with drug AUC values from cell-line experiments. For target prioritization, we retrieved protein targets of drugs from chEMBL, and performed target enrichment analysis on the sRGES-ranked drug list. The enrichment score was calculated using the single sample gene set enrich analysis (ssGSEA), and its significance was computed by a permutation test.

Microenvironment deconvolution

We performed cell type deconvolution of bulk RNA-seq data using both CIBERSORTx [32] (estimating 22 cell types) and xCell [33] (estimating 64 cell types), implemented using the R package immunedeconv [32]. Additionally, we inferred receptor-ligand interaction activity from bulk transcriptomes using the absolute model of BulkSignalR R package [34], followed by univariate Cox regression to evaluate the prognostic significance of each receptor-ligand pair. Corrected p-values were calculated with the BH method.

Single-cell transcriptomic preprocessing

We retrieved a ESCC single-cell RNA sequencing atlas from GEO and processed it into Seurat object using the Seurat R package [35] for downstream analysis. Cell subtype marker genes were retrieved from corresponding paper.

For quality control, we retained cells with ≥200 features and <10% mitochondrial gene, while removing genes expressed in <0.1% of cells along with all mitochondrial and ribosomal genes, followed by doublet detection and removal using the DoubletFinder R package [36]. Count matrices were log-normalized (scale factor = 10,000) via the NormalizeData function, and the top 2,000 highly variable genes (HVGs) were identified using FindVariableFeatures function with the vst method. HVG expression matrices were then scaled and centered using ScaleData. Dimensionality reduction was performed via PCA with optimal PC numbers determined by the ElbowPlot method, followed by batch effect correction using Harmony. These selected PCs were subsequently employed for cell clustering, constructing sNN graphs (FindNeighbors) followed by Louvain algorithm-based clustering (FindClusters) and UMAP visualization.

For cell annotation, we first assigned major cell types based on expression of canonical marker genes. Subsequently, using cell subtype markers provided by Zhang et al., we computed AUC scores for each subtype within major cell categories via the AUCell R package [37], ultimately assigning each cell its optimal subtype based on maximum AUC values.

We quantified gene module activity for each cell using the AddModuleScore function. Epithelial scores were computed as the average expression of canonical epithelial markers (KRT14, KRT17, KRT6A, KRT5, KRT19, KRT8, KRT16, KRT18, KRT6B, KRT15, KRT6C, KRTCAP3, SFN, EPCAM), while mesenchymal scores represented the mean expression of established mesenchymal markers (VIM, CDH2, FOXC2, SNAI1, SNAI2, TWIST1, FN1, ITGB6, MMP2, MMP3, MMP9, SOX10, GSC, ZEB1, ZEB2, TWIST2). The epithelial-mesenchymal transition (EMT) score was derived by subtracting epithelial scores from mesenchymal scores. Proliferation scores were calculated based on the average expression of S-phase and G2M-phase genes from the cc.genes function. Additionally, we assessed the differentiation potential using the CytoTRACE2 R package [38].

Cell trajectory inference

Monocle2 R package [39] was used to perform cell trajectory inference. Monocle2 generates the trajectory using the principal graph algorithm. The top 50 highly expressed genes were selected for cell ordering. Dimensionality reduction and trajectory construction were executed via the reduceDimension function, followed by visualization of gene expression pattern on pseudo-time axis using the plot_genes_in_pseudotime function.

Definition of EMT macro states

Using the mclust R package [40], we modeled predicted pseudo-time scores as a mixture of three Gaussian distributions, with curve intersection points serving as data-driven boundaries to classify EMT progression into three distinct states: EMT-early, EMT-stable, and EMT-late.

Inference of metastatic potential and survival/functional dependency

We employed a K-nearest neighbors (KNN) approach to project cell line transcriptomic data onto single-cell differentiation trajectories. Batch effect correction was performed by merging cell line and single cell count matrices followed by the ComBat normalization from sva R package. Using the get.knn function from FNN R package [41], each cell line was mapped to its k-nearest single-cell neighbors along the trajectory. Through systematic evaluation of k-values, we determined k≥50 reliably produced three distinct peaks in the pseudo-time density distribution. Cell line pseudo-time estimates were derived from median values of the 50 nearest single cells, with EMT phase assignment based on established pseudo-time boundaries, thereby facilitating identification of phase-specific metastatic competencies and survival-associated genes.

We extracted the gene knock down profiles generated by shRNA interference across 4,371 target genes and 20 cell lines from LINCs. However, due to the absence of ESCC cell lines, we employed a connectivity method to generate simulated perturbation profiles. We first generated the disease signature using tumor and adjacent normal samples from our bulk RNA-seq cohorts (XJ and TCGA). Subsequently, we computed connectivity scores between disease signature and each perturbation signature using the connectivity map (CMAP) method. Finally, for perturbation signatures across multiple treatment cell lines, concentrations, and times, the signature closest to the median connectivity score was selected to assign a unique signature for each knock down gene.

Spatial transcriptomic preprocessing

We first loaded the feature-barcode matrix and corresponding high-resolution tissue image using the Load10X_Spatial function. Data normalization was performed via SCTransform with retention of the top 2,000 variable features, followed by PCA dimensionality reduction for downstream analyses. For quality control, the localOutliers function from the SpotSweeper R package [42] was used to identify and remove low quality spots: (i) low library size, (ii) few detected unique features, or (iii) high mitochondrial gene content percentage.

Inference of copy number alterations

To delineate the tumor regions within our spatial transcriptomics dataset, we used the STARCH Python package [43] designed to infer copy number alterations (CNAs). STARCH identifies tumor clones (setting K = 2 clones) and non-tumor spots. It confirms the identification of normal spots by clustering the first principal component into two clusters using K-means. Changing the value of K alters the number of identified tumor clones, but the number of cells labeled as tumor cells remains the same. We then annotated tissue regions in conjunction with manual annotations from the pathologist (Z.N. and Q.Z.).

Spot deconvolution

We performed cellular deconvolution of spatial transcriptomics dataset using the spacexr R package [44], which requires both spatial transcriptomic data (including gene expression counts and spatial coordinates) and reference single-cell RNA sequencing data with cell type annotations as input. We ran spacexr in doublet mode to predict one or two predominant cell types for each spot while simultaneously generating a weight matrix representing the proportions of all possible cell types. We first conducted cell subtype deconvolution across all spots, then performing EMT-state-specific deconvolution focused on spots containing mesenchymal-like epithelial cells.

Spatial colocalization analysis

Using a KNN approach, we first identified neighboring spots surrounding each interested spots, then expanded concentric circular zones (radius levels 1-5) outward from the center to define spatially graduated neighborhoods. Comparative analysis of gene/pathway expressions across these zonal partitions revealed distance-dependent patterns: increased expression indicated spatial co-localization with the interested spots, while decreased expression suggested spatial exclusion.

Development and validation of the deep learning model

To build sample labels, we categorized 175 specimens into high, medium, and low CACNA1C expression groups through Gaussian mixture modeling. The corresponding H&E-stained whole-slide images were randomly divided into training and test sets while maintaining a 4:1 ratio, with reproducibility ensured by setting the random seed to 42. Model evaluation employed five-fold cross-validation, where performance metrics were averaged across all test folds.

For digital pathology image management, we established a standardized preprocessing workflow beginning with background correction using combined Gaussian filtering and OTSU, followed by tissue segmentation into non-overlapping 256×256 patches at 40×magnification. Image normalization was subsequently performed using Reinhard method with Z-score standardization of RGB channels to ensure intensity consistency.

Our analytical framework incorporated four attention-based multiple instance learning (AMIL) architectures [45] (Attention_mil, Clam_sb, Clam_mb, TransMIL), with feature extraction initialized using CTransPath pretrained weights [46]. The models were trained for 32 epochs using cross-entropy loss function and Adam optimizer with an initial learning rate of 1×10⁻⁴, with the best-performing model from cross-validation selected for subsequent analysis.

To investigate the clinical relevance of model attention patterns, we extracted patch-level features from the optimal AMIL model's average pooling layer and conducted unsupervised clustering analysis using the Seurat package. Computational efficiency was maintained by subsampling 200 representative patches per specimen for dimensional reduction and cluster visualization.

Results

Pathological image analysis reveals an EMT-driven prognostic subtype

Clinically relevant cell types, distributions, and interactions in pathological images were characterized by deep learning-based cell segmentation and spatial mapping. Following strict inclusion criteria (see Methods), we retrospectively obtained 254 H&E-stained whole slide images (WSIs) of primary surgical resection specimens from 254 ESCC patients across three cohorts (XJ = 125, GDPH = 78, TCGA = 50) along with their clinical information (Figure 1A). The WSIs were segmented into 4096×4096 patches at 40x magnification. For each patch, cell segmentation was performed to identify tumor cells (T), immune cells (I), and other stromal cells (S), and a local cell interaction map was constructed using the minimum spanning tree (Figure 2A). Edges longer than 50 pixels were discarded, and up to five edges were constructed per cell (denoted from v1 to v5). We quantified the distribution of edge length using six statistical methods (minimal, maximum, mean, variance, skewness, and kurtosis), generating a multiparameter phenotypic descriptor for each patch. These patch-level descriptors were further aggregated into WSI-level features using four statistical measures (mean, variance, skewness, and kurtosis) for subsequent sample-level analysis.

To investigate the clinical relevance of pathological features, we assessed their associations with overall survival (OS) in the XJ cohort via a univariate Cox regression model (Figure 2B, Table S2). Five types of cell interaction features (I-I, I-S, S-S, T-I, and T-S) were proven to be independent prognostic factors for ESCC. By combining feature interaction and aggregation types, we performed enrichment analysis and noted that immune cell features (I-I_sd, I-S_kurtosis, and I-S_sd) correlated with favorable survival, whereas stromal cell features (S-S_sd, T-S_kurtosis) correlated with poorer outcomes (Figure S1A-B) and had the greatest variability across samples (S-S_skewness, T-S_skewness) (Figure S1C, Table S3). Correlation analysis revealed strong associations between the same feature type and weak associations across different types (Figure 2C), and the T-S features demonstrated the strongest independence. Previous studies have shown that the interplay between cancer-associated fibroblasts (CAFs) and epithelial cells can promote malignant transformation and contribute to the formation of an immunosuppressive TME in ESCC [47, 48]. The 65 significant prognostic features were then used to construct a LASSO-Cox regression proportional hazards model to obtain a continuous risk score. This model significantly stratified patients into two risk groups: the XJ cohort (HR = 2.75, 95% CI = [1.70, 4.45], log-rank p < 0.0001), the TCGA cohort (HR = 2.55, 95% CI = [1.07, 6.10], log-rank p = 0.022), and the GDPH cohort (HR = 2.21, 95% CI = [1.05, 4.65], log-rank p = 0.039) (Figure S1D-F). To analyze the biological basis of stratification, we calculated the correlation between the risk scores and the gene set variance analysis (GSVA) scores of the hallmark collection. The top 5 pathways were the coagulation cascade, KRAS upregulation, apoptosis, angiogenesis, and epithelial-mesenchymal transition (Figure S1G), which play important roles in the development of ESCC [49-53].

To further explore the associations of pathological features with clinical prognosis and biological molecules, we merged the XJ and TCGA cohorts (n = 175) with the corresponding bulk-seq profiles. After unsupervised clustering (see Methods), we identified two pathological subtypes with significant overall survival differences (HR = 2.66, 95% CI = [1.75, 4.044], log-rank p = 0.0065) (Figure 2D-E). Mirroring our findings from the univariate Cox regression, subtype 1 exhibited elevated stromal cell features and short survival times, whereas subtype 2 presented increased immune cell features and prolonged survival (Figure 2F-G). Differential expression analysis from bulk RNA-seq revealed 754 significantly overexpressed genes in subtype 1 and 1,035 in subtype 2 (|Log2FC| ≥ 1, p ≤ 0.05) (Figure 3A). Overrepresentation analysis (ORA) of Gene Ontology terms revealed that cell adhesion, collagen formation, epithelial cell proliferation, and differentiation were specifically upregulated in subtype 1, whereas leukocyte immunity, T-cell differentiation and regulation were upregulated in subtype 2 (Figure S2A). Moreover, gene set enrichment analysis (GSEA) of hallmark samples revealed that subtype 1 was enriched in the EMT and hypoxia pathways in addition to the common P53 and TGF-β pathways, and subtype 2 was enriched in the immune response pathway (Figure 3B). Deconvolution analysis also revealed that the fractions of keratinocytes and immunosuppressive M2 macrophages were increased in subtype 1, whereas those of B cells and T cells were increased in subtype 2 (Figure S2E-F). Whole-exome sequencing (WES) analysis revealed 5 significant somatic mutations (Fisher's exact test p < 0.05) (Figure S2C), including TTN (p = 0.033, 95% CI = [1.00, 23.22]), PKHD1L1 (p = 0.020, 95% CI = [0.04, 0.94]), EHBP1 (p = 0.015, 95% CI = [0.001, 0.81]), LRRC56 (p = 0.015, 95% CI = [0.06, 0.81]), and ZNF429 (p = 0.0043, 95% CI = [0, 0.47]). However, the tumor mutation burden (TMB) and fraction of genome altered (FGA) metrics were not significantly different between the subtypes (Figure S2D).

 Figure 1 

Schematic workflow of this study. (A) The hematoxylin & eosin (H&E)-stained whole-section images (WSIs) of esophageal squamous cell carcinoma (ESCC) used in this study are part of routine clinical practice. We used a cell segmentation model to extract pathohistological features from each WSI. Unsupervised clustering revealed pathological subtypes with prognostic value. Bulk-seq data elucidated the molecular dynamics and tumor microenvironment (TME) landscapes of each subtype, allowing us to screen functional ligand-receptor pairs and therapeutic targets in combination with publicly available data. EMT, epithelial-to-mesenchymal transition. (B) Reconstruction and validation of ESCC-specific EMT trajectories. Using single-cell pseudotime analysis, we delineated three macro-EMT states. We combined cell line experimental data to identify core molecules that play central roles in EMT progression and validated their spatial colocalization via spatial transcriptome data. (C) Building a deep learning pathology model for predicting the expression and distribution of molecular targets. We elucidated its interpretability by unsupervised clustering of and matching the predictive attention heatmaps with in situ spatial transcriptome data.

 Figure 2 

Pathological analysis reveals prognostic phenotypes in esophageal squamous cell carcinoma (ESCC). (A) Local nuclear graph visualization of the tumor cells, immune cells, stroma cells and edges among them. (B) Significant univariable associations of pathological features with overall survival in the XJ cohort. (C) Heatmap of correlations between significant prognostic features in the XJ cohort. (D) Workflow of unsupervised clustering. (E) Principal component visualization of the unsupervised clustering result from the merged dataset. (F) Kaplan-Meier analysis of the pathological subtypes. (G) Top 50 contributing pathological features in the principal component analysis (PCA).

 Figure 3 

Tumor microenvironment (TME) and therapeutic characteristics in ESCC prognostic risk stratification. (A) Volcano plot of differentially expressed genes between pathological subtypes, highlighting the top 20 upregulated and downregulated genes. (B) GSEA of hallmark gene sets. (C) Spearman correlation between the drug response AUC and sRGES in the high-risk subtype representative cell line EFO27. (D) Core receptor-ligand pairs. (E) Heatmap of correlation z score between the predicted pathway scores and GSVA scores. (F-G) UMAP visualization of major cell types and subtypes identified in the Zhang_2021 dataset. (H) Heatmap depicting the cellular subtype localization of core receptor-ligand pairs. (I) Expression heatmap of MES-specific receptors and their corresponding ligands. (J) Expression of MES-specific targets across epithelial cell subtypes.

We hypothesized that EMT primarily contributes to the stratification of pathological subtypes. To demonstrate this, we utilized a regularized generalized linear model to learn the associations of pathological features with three types of molecular pathways summarized manually by Jiang et al. [54]. This analysis revealed that the G2M checkpoint pathway in the proliferation type and the EMT pathway in the migration/immune type were more consistent between cohorts than were pathways in the metabolism type (Figure 3C) (Table S4). Overall, these results indicated that epithelial cells exhibiting high proliferation and differentiation potential, an immunosuppressive TME and enrichment of EMT signatures were characteristic of subtype 1, whereas subtype 2 was related to an immune-activated phenotype and favorable clinical outcomes, and EMT primarily drove the prognostic stratification of pathological subtypes.

Vulnerability and dependency across pathological subtypes

For the EMT-enriched & high-risk subtype 1, we next sought to identify potential therapeutic agents and targets using a computational approach. The disease signature was defined as a list of differentially expressed genes (DEGs) between tumor and peritumoral samples of subtype 1. By combining drug perturbation transcriptome profiles of cancer cell lines [31], we screened drugs that can inhibit upregulated genes, stimulate downregulated genes, and eventually reverse the gene expression pattern of the disease signature. We acquired a ranked list consisting of 12,443 small-molecule compounds (Table S5). The VEGF-receptor-2-kinase-inhibitor-IV hit was the top-ranked compound with a known mechanism of action [55, 56]. In addition, we used the transcriptome profiles to select a cancer cell line most related to subtype 1 and then leveraged the published drug sensitivity data of the representative cell line for in silico validation (see Methods). The summarized reverse gene expression scores (sRGES) were significantly positively correlated with the drug efficiency AUCs (Table S6), which increased the confidence of the ranked drug list (Figure 3D). Therefore, we performed enrichment analysis to identify targets whose corresponding drugs were significantly enriched at the top of the prediction. In total, 25 genes were identified as potential targets (NES > 0, p ≤ 0.05) (Table S7). Among these genes, the high expression of CACNA1C, which encodes an alpha-1 subunit of a voltage-dependent calcium channel, has been proven to be associated with poor differentiation of ESCC [57]. Moreover, high expression of its corresponding long-chain noncoding RNA CACNA1C-AS2 inhibits the proliferation, migration, and invasion of esophageal cancer [58]. These findings highlight CACNA1C as a potential therapeutic target in ESCC subtype 1.

Considering that the pathological features were derived from a cell-cell interplay map, we next investigated ligand-receptor interactions at the molecular level. On the basis of the manually curated ligand-receptor interaction database [59], we identified 824 pairs in our dataset, with 381 upregulated in subtype 1 and 441 upregulated in subtype 2 (Table S8). Univariate Cox regression revealed a set of 109 pairs significantly associated with OS, with 97 classified as unfavorable prognostic factors and 12 as favorable prognostic factors. Opposite trends for the hazard ratios (negative vs. positive values) were observed for epidermal growth factor receptors (EGFRs) and fibroblast growth factor receptors (FGFRs), indicating their antagonistic functions (Figure S2G). A summary of the pair-associated pathways revealed that the upregulated pairs in subtype 1 were involved mainly in collagen formation and planar cell polarity, which are also independent unfavorable factors (Figure 2H-I). This analysis enabled us to identify a set of 85 pairs (defined as core pairs, including 42 ligands and 32 receptors) that were positively correlated with the LASSO-Cox risk score, upregulated in subtype 1, and unfavorable factors in terms of prognosis (Figure 3E).

The bulk-seq data includes data on transcriptional programs from a variety of cell types, which could introduce noise signals to our analysis. Therefore, we next validated candidate therapeutic targets and functional ligand-receptor pairs at single-cell resolution (Figure 1B). We incorporated an ESCC single-cell atlas from the Zhang_2021 dataset [60], including 208,659 cells from 64 donors with a total of 128 samples (64 tumors and 64 peritumoral samples). This atlas is represented by 8 major cell types and 51 cell subtypes (Figure 3F-G) (Figure S2J). At the major cell type level, the expression of core pairs was universally upregulated in epithelial cells (Table S8), highlighting their close interactions (Figure S2K). When the signals were refined to the cell subtype level, both the incoming and outgoing signals were dominated by mesenchyme-like epithelial cells (MESs) (Figure 3H) (Figure S2L). Among the core pairs pertaining to epithelial cells (Figure 3I), LAMA3/LAMB3/LAMC2-ITGA3/ITGB4 were significantly enriched in MESs, indicating their autocrine regulation. In addition, we detected high-frequency interactions between CAFs and MESs: CAF4 and CAF2 could interact with ITGA3⁺MES through FN1 and PLAU, respectively, and LAMB1⁺CAF4 could interact with both ITGA3⁺MES and ITGB4⁺MES. Although many genes were not detected, we confirmed that 14 therapeutic candidates were specifically upregulated in epithelial cells (Table S9). Strikingly, CACNA1C was the only significantly highly expressed gene in the MESs (logFC = 2.087, p < 0.001) (Figure 3J). Overall, we identified CACNA1C as a promising therapeutic target for the EMT-enriched subtype 1 subtype and revealed widespread autocrine interactions and paracrine crosstalk of MESs.

Reconstruction and exploration of EMT trajectories

We focused on 5,986 MESs from the Zhang_2021 dataset for further analysis. Unsupervised clustering resulted in seven clusters (Figure 4A). Both the epithelial and mesenchymal programs exhibited cluster-specific enrichment (Figure 4B), resulting in a gradient distribution of EMT scores across clusters (Figure 4C). We inferred the development of MESs by computing a diffusion map and ordering them along a pseudotime axis. Using cluster 0 (with the highest epithelial program) as the origin, trajectory inference predicted two main branches: branch 1 developed through cluster 2 to cluster 1, and branch 2 moved through cluster 3 instead, ending in cluster 4 and cluster 5 (Figure 4D). These trajectories demonstrated that EMT did not follow monotonic linear progression but rather exhibited a stable intermediate state.

To characterize the dominant states governing the continuum of transcriptional activity, we segmented the trajectories on the basis of pseudotime values using a Gaussian mixture model and identified three macro states (Figure 4E): the early state (EMT-early), intermediate/residency state (EMT-stable), and late state (EMT-late). These states were robust to varying pseudotime values, with EMT scores progressively increasing across the three states (Figure 4F) (Figure S3A). Branch-dependent gene analysis revealed regulators governing state transitions (Figure 4G). As expected, keratins (KRT23/81) and collagens (COL1A1/2, COL3A1, COL5A/2, and COL6A3) predominantly shaped epithelial and mesenchymal phenotypes during early and late states, respectively, and TGFBI served as the phenotypic stability factor (PSF) for maintaining the EMT-stable state. We noted that chemokines (CCL21, CXCL8/14) and human leukocyte antigens (HLA-DRB1) were highly expressed across the three states, indicating immune system involvement in EMT. Notably, cells in the EMT-early state exhibited significantly greater proliferative ability, and cells in the EMT-stable state exhibited significantly greater stemness. In terms of their differentiation potential, most MESs maintained lineage-restricted limited differentiation potential (Figure 4H). Recent studies have demonstrated that cells with hybrid E/M phenotypes can either be permanently “locked” in one state or dynamically switch states, which is called epithelial-mesenchymal plasticity (EMP) [61, 62]. These results suggested that the ESCC-specific EMT trajectories not only recapitulated classical marker genes but also captured widespread phenotypic plasticity and crosstalk with the immune system, supporting the development of a reliable, generalizable EMT model.

To investigate how EMT states influence metastatic potential, after removing nonbiological batch effects between in vitro cell lines and in vivo cell models (Figure S3B-C), we projected the corrected metastatic-annotated ESCC cell line profiles onto EMT trajectories (see Methods). Since only epithelial cells were included in the following analysis, confounding effects from mesenchymal cells were precluded. Our analysis robustly captured three EMT macrostates in the cell lines (Figure S3D), with a strong positive correlation between pseudotime values and EMT scores (Figure S3E). As expected, the nonmetastatic group presented relatively lower EMT scores, but the weakly metastatic group presented significantly higher EMT scores than did the metastatic group did (Figure 4H), and we observed a strong negative correlation between EMT scores and metastatic potential across all ESCC cell lines (Figure S3F). This pattern was corroborated by single-cell data, which revealed a decrease in metastatic potential with increasing EMT (Figure 4H) (Figure S3G). The reversal of EMT, a process called mesenchymal-epithelial transition (MET), has been proven to promote metastatic outgrowth at distant sites [63, 64]. These findings underscore the importance of homeostasis between epithelial and mesenchymal programs and reveal that cells in the EMT-early state, which exhibit a predominant epithelial program, have enhanced metastatic potential.

ITGA3 and ITGB4 exhibit functional antagonism in early EMT

We identified distinct expression patterns of MES-specific receptors associated with different EMT states: IGF1R and PLXNA1 were preferentially expressed during the EMT-stable state, whereas ITGA3 and ITGB4 were enriched in the EMT-early state (Figure S3H). Intriguingly, ITGA3 depletion concurrently suppressed epithelial programs and activated mesenchymal programs, whereas ITGB4 depletion induced the opposite effects (Figure 4J). This opposing effect was further evidenced by strongly negatively correlated transcriptome profiles upon perturbation (Figure S3I), indicating fundamental differences in the downstream regulatory consequences. In contrast, IGF1R and PLXNA1 knockout did not significantly affect transcriptome profiles (Figure S4J). Collectively, these findings establish an antagonistic relationship between ITGA3 and ITGB4 in the EMT-early state, where ITGB4 promotes EMT and ITGA3 drives MET. The dynamic balance of these processes determines the differentiation fate of cells in the EMT-early state (Figure 4L).

 Figure 4 

Reconstruction and exploration of epithelial-to-mesenchymal transition (EMT) trajectories in ESCC. (A) UMAP visualization of MES subclusters identified from the Zhang_2021 dataset. (B-C) Expression patterns of phenotype marker genes across MES subclusters. (D) Differentiation trajectories of MESs. (E) Gaussian mixture modeling defining three EMT states (left), and relationships between MES subclusters, pseudotime states, and EMT states (right). (F) Landscape of EMT trajectories. (G) Heatmap of branch-dependent genes during state transitions. (H) Comparison of proliferative capacity (upper-left) and stemness (upper-right) across EMT states, EMT scores across metastatic groups in ESCC cell lines (bottom-left), and metastatic potential across EMT states (bottom-right) via the Wilcoxon signed-rank test. (I) Expression localization (upper) and survival dependency (bottom) of CACNA1C during EMT. A lower dependency score means that a gene is more likely to be dependent on a given cell line. (J) Effects of gene knockdown on the expression of keratins and collagens. (K) GSEA of genes differentially expressed after CACNA1C knockdown using the Reactome gene sets, highlighting the pathways associated with the most significantly up- and downregulated genes. (L) Schematic diagram of CACNA1C-targeting mechanisms.

Inhibition of CACN1C reprograms the malignant phenotype

The phenotypic plasticity of MESs and the dynamic nature of epithelial and mesenchymal programs during metastasis highlight a possible weakness of cell transdifferentiation. As CACN1C was identified as a promising therapeutic target for the poorly prognostic subtype, we next investigated its function in the EMT trajectory. CACNA1C⁺MESs demonstrated a lower expression level in cells in the EMT-early state than in those in the EMT-stable/late states (Figure 4I). Consistently, projection of corrected dependency-annotated ESCC cell line onto EMT trajectories revealed that cells in the EMT-stable/late states showed significantly greater survival dependency than those in the EMT-early state (Figure 4I). We postulated that the elevated expression of CACN1C reflected enrichment of cells in the EMT-stable/late state, potentially explaining its correlation with the poor differentiation of ESCC. To delineate its mechanism, we further analyzed the gene knockdown transcriptome profiles of the ESCC cell lines. We found that CACNA1C depletion concurrently upregulated epithelial programs, downregulated mesenchymal programs (Figure 4J), and significantly suppressed the cell cycle pathway (Figure 4K). Therefore, we hypothesized that targeting CACN1C simultaneously inhibits both the mesenchymal program and the proliferative program in cells in the EMT-stable/late state and initiates their transdifferentiation into normal epithelium-like cells (Figure 4L).

Screening of human ESCC cell line from DepMap illustrated the highest basic CACNA1C expression and dependency in TE-8, which had weakly metastatic potential (Figure 5A). In addition, TE-8 showed sensitivity to more than half of CACNA1C targeted compounds (Cell viability < 0) (Figure 5B). To validate the impact of CACNA1C on the malignant phenotype, we performed in vitro gene knock-down and functional assays using TE-8 (Figure 5C). The mRNA expression of cancer stemness related genes, including BMI1, SOX2, KLF4, MYC, OCT4, and NANOG were analyzed in paired cell line models with or without knock-down. RT-qPCR showed that BMI1 and SOX2 were significantly decrease, and KLF4 was lightly increased compared with controls (Figure 5D). The activation of BMI1 and SOX2 has been reported to be drivers of ESCC [65]. KLF4 is enriched in normal esophageal epithelium, and it has been shown that its expression is lost in ESCC and associated with poor prognosis [66], thus the upregulated of KLF4 might indicate the tendency of recovery of cell phenotype from malignant to normal. In contrast, the expression of other stemness genes, including MYC, OCT4, and NANOG, remained largely unchanged, which indicated a selective disruption of a specific regulatory module of cancer stemness. Additionally, the chamber invasion assay (Figure 5E) and the colony formation/proliferation assay (Figure 5F) showed that the CACNA1C enhanced the invasion and proliferation process of TE-8. These findings collectively highlight the essential role of inhibiting CACNA1C in reprogramming ESCC malignant phenotype.

Spatial colocalization of EMT niches

Given the extensive interactions among MESs, we hypothesized that their spatial colocalization might establishe EMT niches. We employed high-resolution spatial transcriptomic technologies including Visium HD for discovery (tumor E0: 152,517 bins at 16 μm resolution) and Stereo-seq for verification (tumors E1 and E2: 130,927 and 106,236 bins at 20 bin resolution, respectively). Following a uniform quality control pipeline (Figure S4A) (Figure S5A), we performed spot deconvolution based on single-cell data, which generated a spatial map consistent with the tissue morphology (Figure 6A). Copy number variation (CNV) inference delineated the tumor region, tumor-normal epithelium interface, and normal/stromal regions. With reference to the pathologist-annotated tissue masks, we further defined the tumor-stroma interface as a 3-spot-wide zone adjacent to the tumor boundary.

Oxidative phosphorylation-characterized epithelial cells (OXDs) and MESs were the predominant cell types in the spatial map, and they exhibited distinct zonation patterns. OXDs were localized in the tumor nest, whereas MESs surrounded them in the periphery (Figure 6A). The transcriptomic similarity between MESs and OXDs was corroborated by single-cell UMAP projection (Figure 3I). The application of EMT trajectories to spots containing MESs revealed a reliable spatial map of EMT states (Figure 6B). Intriguingly, we noted the spatial overlap of EMT-stable cells and OXDs. Reactive oxygen species (ROS) accumulate during EMT initiation and cancer progression [67], and this spatial association suggests that OXDs may represent an alternative origin distinct from that of EMT-early cells. A quantitative comparison of the cell proportions in tissue regions demonstrated the enrichment of cells in the EMT-early state at tumor-normal interfaces, while cells in the EMT-late state predominated at tumor-stroma interfaces, emphasizing the stroma-oriented invasion of cells in the EMT-late state (Figure 6C).

 Figure 5 

Inhibition of CACN1C represses stemness, invasion and proliferation of TE-8 human ESCC cell line. (A) Basic characteristics of ESCC cell line. (B) ESCC cell line viability after disturbance of CACNA1C targeted compounds. (C) Scheme of cell line functional experiment. (D) The relative expression of CACNA1C and stemness markers in TE-8 cell line expressing siRNAs targeting CACNA1C (siRNA-CACNA1C) or control siRNA (siRNA-NC). Quantitative PCR with reverse transcription (RT-qPCR) data are normalized to GAPDH mRNA expression. n = 6 replicates per group. Statistical analysis was performed using two-sided Student's t-test. (E) Left, Crystal violet staining imaging of invasion assay. Right, Quantification showing the knock-down effect on TE-8 invasion. (F) Left, Crystal violet staining imaging of colony formation assay. Right, Quantification showing the knock-down effect on TE-8 clone formation and proliferation.

 Figure 6 

Spatial cell colocalization along the EMT trajectories. (A) Spatial mapping of a VisiumHD sample (E0) with 16-μm bins colored by cell type predicted by deconvolution via the single-cell reference dataset (right) and MES-containing spots colored according to tissue region (bottom). (B) Spatial mapping of EMT states. (C) Stacked histogram with the proportions of EMT states across tissue regions. (D) Spatial visualization of MES-specific receptors. (E) Spatial visualization of ITGA3 and ITGB4 in MES-containing spots. (F) Spatial correlation between ITGA3/ITGB4-specific EMT-early spots and regional EMT activity. (G) Schematic of workflows for Xenium experiment. (H) Top, Diagram showing the quantification of spatial proximity of cells, and violin and bar plot showing EMT-early subtype classification. Bottom, Kernel density plots displaying the EMT-early distribution along X and Y axis, and stacked bar plots showing the proportion of EMT-early. P values calculated by Mann-Kendall trend test.

Spatial mapping of MES-specific receptors showed their predominant localization within the tumor region (Figure 6D). Specifically, ITGA3 exhibited preferential expression toward the normal epithelium, whereas ITGB4 presented tumor- and stroma-oriented expression (Figure 6E). Spots containing EMT-early state cells were classified into four groups based on ITGA3 and ITGB4 expression (Figure S4B). As predicted, the ITGA3⁺ITGB4^- spots had the lowest EMT scores, whereas the ITGA3^-ITGB4⁺ spots showed significantly higher scores (Figure 6F). Despite the limited spots detected (n = 407), ITGA3⁺ITGB4⁺ spots presented the highest EMT scores, which was recapitulated by the single-cell data (Figure S4C). These findings collectively suggest a cooperative effect in which ITGA3 and ITGB4 are coexpressed in the EMT-early state.

Using a neighborhood enrichment approach (see Methods), we evaluated the environment around ITGA3⁺ITGB4^- and ITGA3^-ITGB4⁺ EMT-early spots (Figure S4E-F). EMT scores decreased progressively with distance from the ITGA3^-ITGB4⁺ spots, but showed the opposite trend for the ITGA3⁺ITGB4^- spots (Figure 6F). Using Stereo-seq data, we consistently observed spatial exclusion of ITGA3 and ITGB4 (Figure S5D-H), their potential cooperative effects, and distance-dependent EMT distributions (Figure S5I-J). In addition, EMT-early specific ligands (LAMA3, LAMB3, LAMC2) showed decreasing expression with distance from both ITGA3⁺ITGB4^- and ITGA3^-ITGB4⁺ spots (Figure S4F-G). In contrast, expression of the CAF4-specific ligand FN1 increased with distance from ITGA3⁺ITGB4^- spots (Figure S4H).

We further employed 10x Xenium technology for in-situ spatial profiling 5,000 mRNAs across 385,959 cells from 6 primary ESCC samples (Figure 6G). After obtaining 14 major cell types using marker genes, we employed label transfer approach to acquire subtypes of epithelial and fibroblast cells using single-cell data as reference (Table S6). Subsequently, we distinguished EMT-early subtypes by the expression levels of ITGA3 and ITGB4, and calculated spatial distances between cells at single-cell resolution. In brief, we not only verified the spatial exclusion between ITGA3⁺ITGB4^- and ITGA3^-ITGB4⁺ EMT-early (Figure 6H), but also the spatial co-localization between ITGA3⁺ITGB4^- EMT-early with CAF4 (Figure S6). Taken together, these results collectively suggest that cells in the EMT-early state colocalize with self-derived or CAF4-secreted ligands to establish EMT niches.

Prediction of CACN1C expression and distribution using pathological H&E images

Deep learning has been successfully applied in various cancer types for extracting clinically relevant features from routine histopathological slides [68-70]. Building on the above findings, we subsequently investigated whether deep learning could accurately predict molecular target expression and whether computational predictions would reveal meaningful clinical and molecular associations.

We initially discretized the CACNA1C expression values into three groups (high/medium/low) as sample labels (Figure 7A), with significant differences in OS between the medium- and low-expression groups, whereas the high-expression group presented a fluctuating Kaplan-Meier curve (Figure 7B), highlighting the biological complexity of ESCC prognosis. We subsequently trained several attention-based multiple instance learning (AMIL) models for three-class classification (Figure 7C) (Figure S7), among which Clam_mb demonstrated superior predictive performance in our dataset, achieving a mean fivefold cross-validated area under the receiver operating characteristic curve (AUROC) of 0.668, outperforming other AMIL models (Clam_sb = 0.658, Attention_mil = 0.648, Trans_mil = 0.628).

To decipher the clinical implications of model predictions, we extracted patch-level features from the average pooling layer of the Clam_mb model, reduced feature dimensionality to 50 using PCA, and identified six patch clusters using Louvain clustering method (Figure 7D). These clusters demonstrated differential attention score contributions: cluster 5 showed the highest attribution for low-expression classification, cluster 0 for medium-expression classification, and cluster 1 for high-expression classification. Whole-slide visualization and H&E-correlated annotation revealed distinct histopathological identities: cluster 0 represented stromal tissue, cluster 1 represented the tumor region, cluster 2 represented muscle/normal epithelium, cluster 3 represented the tumor stroma, and clusters 4 and 5 represented invasive tumor margins (Figure 7E). Collectively, these computationally derived clusters presented unique histopathological signatures, suggesting the tissue-specific expression of CACNA1C in ESCC.

 Figure 7 

Development and validation of the pathological deep learning model. (A) Sample classification on the basis of CACNA1C expression. (B) Kaplan-Meier analysis of the sample classification. (C) Receiver operating characteristic (ROC) curves showing the performance of the Clam_mb model in fivefold cross-validation. (D) UMAP visualization of patch-level clusters derived from deep learning features, with box plots quantifying each cluster's contribution to attention scores. (E) Spatial distribution mapping of computational patch clusters across whole-slide images. (F) Spatial concordance between the model attention score and in situ CACNA1C expression in tumor tissues.

To gain further insights into the in situ relationships between the model predictions and the biological ground truth, we performed Visium sequencing on two ESCC FFPE tissue sections (capturing 4,820 and 4,468 spots from tumors E3 and E4, respectively), with tumor/stroma demarcation through CNV inference. Despite the inherent tissue morphology variations between serial sections, application of the Clam_mb model to adjacent H&E-stained sections revealed strong concordance between attention heatmaps and spatial transcriptomics, where model-identified high-attention regions exhibited corresponding gene-specific upregulation (Figure 7F). Notably, we did not incorporate tissue masks or region-specific cropping during training, and the model selectively focused on high-expression regions in tumor areas while filtering out stromal signals (Figure S7D), demonstrating that our deep learning framework can sensitively detect in situ CACNA1C expression patterns in tumor tissue from H&E images.

Discussion

Understanding the mechanisms by which cell crosstalk drives cancer aggression and metastasis is critical for improving the diagnostic precision and development of targeted therapies for ESCC. Here, we described the development, validation, and explanation of pathological cell-cell interplay features for prognosticating overall survival time using H&E-stained images. We combined a deep learning model for cell segmentation and classification, a minimum spanning tree for topological feature construction and extraction, a Cox regression model and unsupervised clustering to maximize interpretability. We developed and validated our features using three independent cohorts. Our results revealed that image-derived cell interaction features serve as strong predictors of survival outcomes in ESCC patients and are concordant with known molecular pathways. The clinical and biological relevance of these features support their utility and generalizability.

Using combined bulk-seq data, we demonstrated the existence of a high-risk pathological subtype characterized by an EMT-enriched state and emphasized the concomitant immunosuppressive microenvironment. In addition, the ligand-receptor pairs enriched in pathological subtypes demonstrated EMT-promoting effects, with predominant localization between mesenchymal cell-like epithelial cells (MESs) and between MESs and CAFs, indicating both autocrine and paracrine regulation of EMT. Through systematic target screening, we identified CACNA1C as a specifically overexpressed target in MESs, which is consistent with its established association with poor ESCC differentiation.

EMT progression has been proven to be a continuous phenotypic process characterized by multiple intermediate states [71]. However, to our knowledge, relevant investigations have relied on pancancer cell line models, leaving a critical gap in single-cell characterization using clinical specimens. We reconstructed EMT trajectories using the largest available single-cell atlas of clinical ESCC samples. Gradient shifts between epithelial and mesenchymal programs across MESs demonstrated phenotypic continuity, and pseudotime analysis revealed a plateau state during EMT progression. On the basis of pseudotime score distributions, we defined three macrostates for MESs (EMT-early, EMT-stable, and EMT-late), with progressively increasing EMT scores across states. As anticipated, keratins and collagens mainly shape epithelial and mesenchymal phenotypes in the early and late states, respectively, whereas TGFBI functions as a phenotypic stability factor during the stable phase. Notably, chemokines and human leukocyte antigens exhibit stage-dependent expression, indicating that there is crosstalk between EMT and the immune system. Furthermore, cells in the EMT-early state displayed the greatest proliferative capacity, cells in the EMT-stable state presented the highest differentiation potential, and most MESs maintained lineage-restricted limited differentiation potential, collectively demonstrating remarkable phenotypic plasticity during EMT. The full spectrum of EMT intermediate states remains to be fully characterized, and the mechanism of state switches needs to be better understood. Here, the ESCC-specific EMT trajectory delineates three well-supported macrostates reflecting both intrinsic cancer cell alterations and microenvironment adaptations.

It has been reported that cells with hybrid E/M phenotypes generate progeny cells that are either mesenchymal or epithelial and are more prone to migrate [72]. To evaluate the impact of EMT on ESCC metastasis, we projected metastatic ESCC cell lines onto EMT trajectories. This analysis revealed that cells in the EMT-early state had the strongest metastatic potential, highlighting the critical balance between the epithelial and mesenchymal programs. While enrichment of the epithelial program is required for ESCC metastasis, we demonstrated that targeting CACNA1C simultaneously suppresses both the mesenchymal program and cell cycle progression in EMT-stable/late states, driving their transdifferentiation into normal epithelium-like cells. Furthermore, we established that ITGB4 overexpression promotes EMT, whereas ITGA3 enhances MET. Importantly, by exclusively analyzing epithelial cells, we eliminated potential confounding effects from mesenchymal cell contamination.

Recent studies have demonstrated that hybrid E/M cells preferentially localize to the tumor invasive front [73-75]. Similarly, our spatial transcriptomic data revealed enrichment of EMT-early cells at the tumor-normal epithelium interface, whereas EMT-late cells were predominant at the tumor-stroma border. Notably, we observed spatial overlap between EMT-stable and epithelial cells with oxidative phosphorylation characteristics (OXDs), suggesting that OXDs may serve as alternative EMT-initiating cells distinct from those with early EMT. Previous work revealed that ITGA3 and ITGB1 are specifically overexpressed in the tumor-specific keratinocytes (TSKs) of cutaneous squamous cell carcinoma, where they function as receptors for ligands expressed by TAMs and MDSCs [76]. Our findings not only confirmed the mutually exclusive spatial expression of ITGA3 and ITGB4 but also revealed a potential cooperative role in driving EMT progression when they are coexpressed in the EMT-early state.

End-to-end deep learning models have yielded impressive results for diagnostic applications such as the detection of cancer and the prediction of the primary origin of metastases [77, 78]. Attention-based MIL methods learn from patient-level labels and offer explanations in the form of saliency heatmaps that localize relevant regions. Inspired by Calderaro et al. [79], we used spatial transcriptomes as the gold standard for the molecular detection problem. We developed a target prediction model in which computationally derived patch clusters exhibited distinct histopathological signatures, with CACNA1C-high expression clusters specifically localized to tumor regions. Spatial validation demonstrated strong concordance between the model predictions and the in situ biological ground truth. Notably, without incorporating tissue masks or region-specific cropping during training, the model autonomously attended to tumor-specific high-expression zones, confirming both the tissue-specific expression of CACNA1C and the translational potential of this model.

This study has some limitations: (i) As a retrospective analysis of patient cohorts, this study requires future prospective validation to confirm its predictive findings. (ii) Although local cell interaction maps effectively correlate with clinical and biological features, the absence of detailed tissue architecture and refined cellular subtypes may limit analytical depth. This limitation could be addressed through dedicated tissue and cell segmentation models. (iii) ESCC-specific EMT trajectories warrant comprehensive validation in clinical metastatic samples, and the time-course dynamics should be further explored. (iv) Model performance should be improved by several approaches, including expanding training datasets, employing more robust pre-training encoders, obtaining accurate labels via in situ hybridization, combining multiple pathological images from the sample patient and multimodal feature engineering.

Overall, the results of this study systematically delineate EMT dynamics in ESCC through integrative analysis of histopathology, bulk-seq, single-cell transcriptomics, and spatial transcriptomics. By analyzing cellular interaction features from pathological images, we first demonstrated that pathological patterns are driven by EMT dynamics. We further validated the continuum of intermediate EMT states, revealed the pivotal role of the EMT-early state in metastasis and further identified CACNA1C as a therapeutic target for cells in the EMT-stable/late state. These findings provide novel insights into the cellular and molecular mechanisms underlying tumor invasion and metastasis. Finally, we developed an end-to-end deep learning model that predicts therapeutic targets from pathology images, enabling clinically translatable risk stratification and personalized therapy for ESCC patients.

Abbreviations

ESCC: esophageal squamous cell carcinoma; TME: tumor microenvironment; EMT: epithelial-mesenchymal transition; WSI: whole slide image; OS: overall survival; CAFs: cancer-associated fibroblasts; MESs: mesenchyme-like epithelial cells; MET: mesenchymal-epithelial transition; CNV: copy number variation; AMIL: attention-based multiple instance learning.

Supplementary Material

Supplementary figures.

Supplementary tables.

Acknowledgements

We thank Biorender platform for the drawing materials. We thank OE Biotech Co., Ltd (Shanghai, China) and Gene Denovo Honour Biotechnology Co., Ltd (Guangzhou, China) for providing RNA sequencing.

Funding

National Science Foundation for Young Scientists (82302299 and 82202840), Natural Science Foundation for Distinguished Young Scholars of Guangdong Province (2023B1515020043), Top Talent Project of Scientific and Technological Innovation in Tianshan Talents Training Plan of Xinjiang Uygur Autonomous Region (2023TSYCCX0058), Distinguished Young Project of Natural Science Foundation of Xinjiang Uygur Autonomous Region (2022D01E76).

Data and code availability

All public data sources employed are described above. The analysis code can be found in our GitHub repository: https://github.com/AllenSpike/Pathogenomics_EMP.

Author contributions

RZ Chen and CY Xie designed this study, performed analysis, constructed model, and wrote the original manuscript. ZY Ning conducted pathological annotations. M Yang, JH Chen, KH Du, Yihuai Hu, QL Zhang contributed to the multi-center data collection. Zezhuo Su, SJ Zhang reviewed and edited the manuscript. C Han, M Liu, and ZY Liu supervised the study. All authors read and approved the manuscript.

Competing Interests

The authors have declared that no competing interest exists.

References

1. Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I. et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2024;74:229-63

2. NCCN Clinical Practice Guidelines in Oncology (NCCN Guidelines®) for Esophageal and Esophagogastric Junction Cancers. 4th ed: © National Comprehensive Cancer Network, Inc.; 2025

3. Liu J, Xie X, Zhou C, Peng S, Rao D, Fu J. Which factors are associated with actual 5-year survival of oesophageal squamous cell carcinoma? Eur J Cardiothorac Surg. 2012;41:e7-11

4. Vellayappan BA, Soon YY, Ku GY, Leong CN, Lu JJ, Tey JC. Chemoradiotherapy versus chemoradiotherapy plus surgery for esophageal cancer. Cochrane Database Syst Rev. 2017;8:CD010511

5. Yang H, Wang F, Hallemeier CL, Lerut T, Fu J. Oesophageal cancer. Lancet. 2024;404:1991-2005

6. Rhim AD, Mirek ET, Aiello NM, Maitra A, Bailey JM, McAllister F. et al. EMT and dissemination precede pancreatic tumor formation. Cell. 2012;148:349-61

7. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S. et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017;171:1611-24.e24

8. Kalluri R, Weinberg RA. The basics of epithelial-mesenchymal transition. J Clin Invest. 2009;119:1420-8

9. Pastushenko I, Blanpain C. EMT Transition States during Tumor Progression and Metastasis. Trends Cell Biol. 2019;29:212-26

10. Goetz H, Melendez-Alvarez JR, Chen L, Tian XJ. A plausible accelerating function of intermediate states in cancer metastasis. PLoS Comput Biol. 2020;16:e1007682

11. Pastushenko I, Brisebarre A, Sifrim A, Fioramonti M, Revenco T, Boumahdi S. et al. Identification of the tumour transition states occurring during EMT. Nature. 2018;556:463-8

12. Feng X, Shu W, Li M, Li J, Xu J, He M. Pathogenomics for accurate diagnosis, treatment, prognosis of oncology: a cutting edge overview. J Transl Med. 2024;22:131

13. Wang H, Jiang Y, Li B, Cui Y, Li D, Li R. Single-Cell Spatial Analysis of Tumor and Immune Microenvironment on Whole-Slide Image Reveals Hepatocellular Carcinoma Subtypes. Cancers (Basel). 2020 12

14. Yoo SY, Park HE, Kim JH, Wen X, Jeong S, Cho NY. et al. Whole-Slide Image Analysis Reveals Quantitative Landscape of Tumor-Immune Microenvironment in Colorectal Cancers. Clin Cancer Res. 2020;26:870-81

15. Zhao S, Chen DP, Fu T, Yang JC, Ma D, Zhu XZ. et al. Single-cell morphological and topological atlas reveals the ecosystem diversity of human breast cancer. Nat Commun. 2023;14:6796

16. Xu H, Usuyama N, Bagga J, Zhang S, Rao R, Naumann T. et al. A whole-slide foundation model for digital pathology from real-world data. Nature. 2024;630:181-8

17. Xiang J, Wang X, Zhang X, Xi Y, Eweje F, Chen Y. et al. A vision-language foundation model for precision oncology. Nature. 2025;638:769-78

18. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D. et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44:e71

19. Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin AA, Kim S. et al. The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature. 2012;483:603-7

20. Iorio F, Knijnenburg TA, Vis DJ, Bignell GR, Menden MP, Schubert M. et al. A Landscape of Pharmacogenomic Interactions in Cancer. Cell. 2016;166:740-54

21. Jin X, Demere Z, Nair K, Ali A, Ferraro GB, Natoli T. et al. A metastasis map of human cancer cell lines. Nature. 2020;588:331-6

22. Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X. et al. A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell. 2017;171:1437-52 e17

23. Arafeh R, Shibue T, Dempster JM, Hahn WC, Vazquez F. The present and future of the Cancer Dependency Map. Nat Rev Cancer. 2025;25:59-73

24. Graham S, Vu QD, Raza SEA, Azam A, Tsang YW, Kwak JT. et al. Hover-Net: Simultaneous segmentation and classification of nuclei in multi-tissue histology images. Med Image Anal. 2019;58:101563

25. Gamper J, Alemi Koohbanani N, Benes K, Graham S, Jahanifar M, Khurram SA. et al. PanNuke Dataset Extension, Insights and Baselines. 2020. p. arXiv:2003.10778.

26. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8:118-27

27. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550

28. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (N Y). 2021;2:100141

29. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7

30. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33:1-22

31. Zeng B, Glicksberg BS, Newbury P, Chekalin E, Xing J, Liu K. et al. OCTAD: an open workspace for virtually screening therapeutics targeting precise cancer patient groups using gene expression features. Nat Protoc. 2021;16:728-53

32. Sturm G, Finotello F, Petitprez F, Zhang JD, Baumbach J, Fridman WH. et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics. 2019;35:i436-i45

33. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biology. 2017;18:220

34. Villemin JP, Bassaganyas L, Pourquier D, Boissière F, Cabello-Aguilar S, Crapez E. et al. Inferring ligand-receptor cellular networks from bulk and spatial transcriptomic datasets with BulkSignalR. Nucleic Acids Res. 2023;51:4726-44

35. Hao Y, Stuart T, Kowalski MH, Choudhary S, Hoffman P, Hartman A. et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol. 2024;42:293-304

36. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst. 2019;8:329-37.e4

37. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G. et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14:1083-6

38. Kang M, Armenteros JJA, Gulati GS, Gleyzer R, Avagyan S, Brown EL. et al. Mapping single-cell developmental potential in health and disease with interpretable deep learning. bioRxiv. 2024

39. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA. et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14:979-82

40. Scrucca L, Fraley C, Murphy TB, Raftery AE. Model-based clustering, classification, and density estimation using mclust in R: Chapman and Hall/CRC; 2023

41. Beygelzimer A, Kakadet S, Langford J, Arya S, Mount D, Li S. FNN: fast nearest neighbor search algorithms and applications. (No Title). 2010

42. Totty M, Hicks SC, Guo B. SpotSweeper: spatially aware quality control for spatial transcriptomics. Nat Methods. 2025;22:1520-30

43. Elyanow R, Zeira R, Land M, Raphael BJ. STARCH: copy number and clone inference from spatial transcriptomics data. Phys Biol. 2021;18:035001

44. Cable DM, Murray E, Zou LS, Goeva A, Macosko EZ, Chen F. et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol. 2022;40:517-26

45. Dolezal JM, Kochanny S, Dyer E, Ramesh S, Srisuwananukorn A, Sacco M. et al. Slideflow: deep learning for digital histopathology with real-time whole-slide visualization. BMC Bioinformatics. 2024;25:134

46. Wang X, Yang S, Zhang J, Wang M, Zhang J, Yang W. et al. Transformer-based unsupervised contrastive learning for histopathological image classification. Med Image Anal. 2022;81:102559

47. Chen L, Zhu S, Liu T, Zhao X, Xiang T, Hu X. et al. Aberrant epithelial cell interaction promotes esophageal squamous-cell carcinoma development and progression. Signal Transduct Target Ther. 2023;8:453

48. Chang J, Lu J, Liu Q, Xiang T, Zhang S, Yi Y. et al. Single-cell multi-stage spatial evolutional map of esophageal carcinogenesis. Cancer Cell. 2025;43:380-97.e7

49. Wang Q, Cao B, Peng L, Dai W, Jiang Y, Xie T. et al. Development and Validation of a Practical Prognostic Coagulation Index for Patients with Esophageal Squamous Cell Cancer. Ann Surg Oncol. 2021;28:8450-61

50. Yu Y, Zhao Q, Cui H, Song L. CCR7 in esophageal squamous cell carcinoma: an identification from single-cell and bulk transcriptome sequencing. Discov Oncol. 2025;16:183

51. Degenhardt K, Mathew R, Beaudoin B, Bray K, Anderson D, Chen G. et al. Autophagy promotes tumor cell survival and restricts necrosis, inflammation, and tumorigenesis. Cancer Cell. 2006;10:51-64

52. Wang S, Liang Y, Zhang J, Wang W, Hong Y, Sun M. et al. The angiogenic genes predict prognosis and immune characteristics in esophageal squamous cell carcinoma: Evidence from multi-omics and experimental verification. Front Oncol. 2022;12:961634

53. Guo D, Sheng K, Zhang Q, Li P, Sun H, Wang Y. et al. Single-cell transcriptomic analysis reveals the landscape of epithelial-mesenchymal transition molecular heterogeneity in esophageal squamous cell carcinoma. Cancer Lett. 2024;587:216723

54. Jiang Y, Sun A, Zhao Y, Ying W, Sun H, Yang X. et al. Proteomics identifies new therapeutic targets of early-stage hepatocellular carcinoma. Nature. 2019;567:257-61

55. Chen L, Lin G, Chen K, Liang R, Wan F, Zhang C. et al. VEGF promotes migration and invasion by regulating EMT and MMPs in nasopharyngeal carcinoma. J Cancer. 2020;11:7291-301

56. Prasad CB, Singh D, Pandey LK, Pradhan S, Singh S, Narayan G. VEGFa/VEGFR2 autocrine and paracrine signaling promotes cervical carcinogenesis via β-catenin and snail. Int J Biochem Cell Biol. 2022;142:106122

57. Shi ZZ, Shang L, Jiang YY, Shi F, Xu X, Wang MR. et al. Identification of genomic biomarkers associated with the clinicopathological parameters and prognosis of esophageal squamous cell carcinoma. Cancer Biomark. 2015;15:755-61

58. Qin Y, Cao T, Xu Y, Jiang X, Ma J, Yu L. Knockdown of lncRNA CACNA1C-AS2 promotes epithelial-mesenchymal transition and enhances the proliferation, invasion and migration of human esophageal cancer cells. Xi Bao Yu Fen Zi Mian Yi Xue Za Zhi. 2023;39:249-57

59. Cabello-Aguilar S, Alame M, Kon-Sun-Tack F, Fau C, Lacroix M, Colinge J. SingleCellSignalR: inference of intercellular networks from single-cell transcriptomics. Nucleic Acids Res. 2020;48:e55

60. Zhang X, Peng L, Luo Y, Zhang S, Pu Y, Chen Y. et al. Dissecting esophageal squamous-cell carcinoma ecosystem by single-cell transcriptomic analysis. Nat Commun. 2021;12:5291

61. Kröger C, Afeyan A, Mraz J, Eaton EN, Reinhardt F, Khodor YL. et al. Acquisition of a hybrid E/M state is essential for tumorigenicity of basal breast cancer cells. Proc Natl Acad Sci U S A. 2019;116:7353-62

62. Pastushenko I, Mauri F, Song Y, de Cock F, Meeusen B, Swedlund B. et al. Fat1 deletion promotes hybrid EMT state, tumour stemness and metastasis. Nature. 2021;589:448-55

63. Ocaña OH, Córcoles R, Fabra A, Moreno-Bueno G, Acloque H, Vega S. et al. Metastatic colonization requires the repression of the epithelial-mesenchymal transition inducer Prrx1. Cancer Cell. 2012;22:709-24

64. Tsai JH, Donaher JL, Murphy DA, Chau S, Yang J. Spatiotemporal regulation of epithelial-mesenchymal transition is essential for squamous cell carcinoma metastasis. Cancer Cell. 2012;22:725-36

65. An L, Li M, Jia Q. Mechanisms of radiotherapy resistance and radiosensitization strategies for esophageal squamous cell carcinoma. Mol Cancer. 2023;22:140

66. Chen Y, Zhu S, Liu T, Zhang S, Lu J, Fan W. et al. Epithelial cells activate fibroblasts to promote esophageal cancer development. Cancer Cell. 2023;41:903-18.e8

67. Jiang J, Wang K, Chen Y, Chen H, Nice EC, Huang C. Redox regulation in tumor cell epithelial-mesenchymal transition: molecular basis and therapeutic strategy. Signal Transduct Target Ther. 2017;2:17036

68. Wang X, Zhao J, Marostica E, Yuan W, Jin J, Zhang J. et al. A pathology foundation model for cancer diagnosis and prognosis prediction. Nature. 2024;634:970-8

69. Lu MY, Chen B, Williamson DFK, Chen RJ, Liang I, Ding T. et al. A visual-language foundation model for computational pathology. Nat Med. 2024;30:863-74

70. Chen RJ, Ding T, Lu MY, Williamson DFK, Jaume G, Song AH. et al. Towards a general-purpose foundation model for computational pathology. Nat Med. 2024;30:850-62

71. Roche J. The Epithelial-to-Mesenchymal Transition in Cancer. Cancers (Basel). 2018 10

72. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144:646-74

73. Dongre A, Rashidian M, Reinhardt F, Bagnato A, Keckesova Z, Ploegh HL. et al. Epithelial-to-Mesenchymal Transition Contributes to Immunosuppression in Breast Carcinomas. Cancer Res. 2017;77:3982-9

74. Bagati A, Kumar S, Jiang P, Pyrdol J, Zou AE, Godicelj A. et al. Integrin αvβ6-TGFβ-SOX4 Pathway Drives Immune Evasion in Triple-Negative Breast Cancer. Cancer Cell. 2021;39:54-67.e9

75. Aggarwal V, Montoya CA, Donnenberg VS, Sant S. Interplay between tumor microenvironment and partial EMT as the driver of tumor progression. iScience. 2021;24:102113

76. Ji AL, Rubin AJ, Thrane K, Jiang S, Reynolds DL, Meyers RM. et al. Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma. Cell. 2020;182:1661-2

77. Campanella G, Hanna MG, Geneslaw L, Miraflor A, Werneck Krauss Silva V, Busam KJ. et al. Clinical-grade computational pathology using weakly supervised deep learning on whole slide images. Nat Med. 2019;25:1301-9

78. Lu MY, Chen TY, Williamson DFK, Zhao M, Shady M, Lipkova J. et al. AI-based pathology predicts origins for cancers of unknown primary. Nature. 2021;594:106-10

79. Calderaro J, Ghaffari Laleh N, Zeng Q, Maille P, Favre L, Pujals A. et al. Deep learning-based phenotyping reclassifies combined hepatocellular-cholangiocarcinoma. Nat Commun. 2023;14:8290

Author contact

Corresponding authors: Meng Liu, MD, PhD (liuxiaomeng198154com). Zaiyi Liu, MD, PhD (liuzaiyiorg.cn).

Citation styles

©2026 Ivyspring International Publisher. Terms of use