A multi-omics integration of bulk and single-cell transcriptomics identifies and validates a 16-gene epithelial-mesenchymal transition prognostic signature in tongue squamous cell carcinoma

Introduction

Tongue cancer represents a prevalent and highly aggressive form of oral squamous cell carcinoma (OSCC), constituting approximately 40–50% of all oral malignancies. It most frequently arises on the lateral borders and ventral surface of the tongue (1). Despite significant progress in the diagnosis and management of oral cancer over recent decades, the long-term outlook for patients with advanced-stage tongue squamous cell carcinoma (TSCC) remains poor, with a five-year survival rate hovering around 50% (2). Current therapeutic modalities, including surgery, radiotherapy, and chemotherapy, are often associated with high rates of recurrence and the development of treatment resistance. This underscores the pressing clinical requirement for reliable prognostic biomarkers and novel therapeutic targets to improve patient survival (3). The epithelial-mesenchymal transition (EMT) has been established as a pivotal biological mechanism implicated in cancer metastasis, therapeutic resistance, and adverse prognosis (4). Prior investigations have emphasized the significant contribution of EMT-associated genes in the pathogenesis of various malignancies, including oral cancers, thereby suggesting the feasibility of developing an EMT-based prognostic model for TSCC (5).

Recent progress in bioinformatics has enabled the integrative analysis of large-scale genomic datasets, allowing for comprehensive exploration of gene expression profiles and cellular interactions in cancer biology (6). This study adopts a multi-omics framework, combining bulk RNA-sequencing (RNA-seq) data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) with single-cell RNA-seq (scRNA-seq) analyses to identify key prognostic genes for TSCC (7). By examining gene expression patterns, the research aims to develop and validate a robust EMT-related gene (ERG) signature for predicting TSCC patient outcomes, while also investigating the associated biological mechanisms and tumor immune microenvironment characteristics (8).

The oncogenic mechanisms of TSCC are complex and not fully elucidated, involving not only epithelial cells but also intricate interactions between stromal and tumor cells (9). The integration of bulk and single-cell transcriptomic data enhances the identification of prognostic markers by enabling the assessment of cellular heterogeneity and microenvironmental influences on tumor progression (10). Furthermore, the application of advanced statistical modeling techniques, such as least absolute shrinkage and selection operator (LASSO) regression, improves the selection of meaningful biomarkers while reducing the risk of overfitting (11). The objectives of this study are to establish a reliable prognostic tool that stratifies patients based on risk profiles and provides insights into the molecular mechanisms of TSCC and its interactions with the immune system (5).

In summary, this research addresses the urgent need for effective prognostic tools in TSCC by utilizing modern bioinformatics techniques to clarify the role of ERGs in disease progression and patient survival. The identification of these biomarkers may lead to improved therapeutic strategies and personalized treatment plans for patients (12). The findings from this study have the potential to significantly advance our understanding of TSCC and enhance clinical decision-making in oncology (13). We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0435/rc).

Methods

Data collection

Gene expression data (RNA-seq transcriptome data) and corresponding clinical information (including gender, age, pathological stage, survival time, and survival status) for TSCC were obtained from the TCGA database (https://portal.gdc.cancer.gov/) within the head and neck squamous cell carcinoma dataset. This dataset included 13 normal control samples and 130 TSCC samples. To enhance the study’s accuracy, a validation dataset (GSE41613) from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) was selected for model verification. GSE41613 contains 97 samples, all of which are OSCC, with no normal samples. Additionally, scRNA-seq data for oral cancer were obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/, accession number: GSE215403), which includes 12 oral cancer samples. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The overall study design and analytical workflow are summarized in Figure 1.

Figure 1 Flowchart of the study. AUC, area under the curve; DEEGs, EMT-related differentially expressed genes; DEGs, differentially expressed genes; EMT, epithelial-mesenchymal transition; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden; TSCC, tongue squamous cell carcinoma.

Screening and analysis of ERGs

Differential gene expression analysis of TSCC data was performed using the DESeq2 package in R, with thresholds set at |log₂ fold change (FC)| >1.5 and P<0.05, using a t-test for group comparisons. A total of 1,373 differentially expressed genes (DEGs) were identified. Volcano plots and heatmaps were generated using R’s ggplot2 and pheatmap packages (displaying the top 25 upregulated and downregulated genes based on absolute fold change values). Subsequently, 4,525 ERGs were identified by searching the term “EMT” on the GeneCards platform (https://www.genecards.org). The intersection of DEGs and ERGs yielded 438 EMT-related differentially expressed genes (DEEGs), which were visualized using a Venn diagram. The top 10 DEEGs with the largest absolute fold changes were selected for group comparison plots.

Identification of ERGs

The TCGA-TSCC dataset served as the training cohort, while a GEO dataset was utilized for validation. Initially, a univariate Cox regression model was implemented via the R survival package to identify 27 genes with statistically significant associations with survival, and their hazard ratios were computed. Subsequently, LASSO-Cox regression analysis was conducted using the glmnet and survival packages, selecting 16 survival-associated genes based on the optimal λ value. This set of genes constituted the ERG prognostic signature for TSCC patients. The R package maftools was employed to generate waterfall plots, depicting the mutation frequency and types of these 16 differentially expressed prognostic genes across samples. Chromosomal localization of these genes was clarified using the RCircos package, and their interrelationships were assessed through correlation analysis performed with the ggcorrplot package.

Development and validation of the prognostic EMT gene signature

Patients with missing survival information were excluded from the analysis to ensure data completeness and reliability of survival outcomes. The genes identified in the previous step were incorporated into a multivariate Cox regression model to further refine the candidate gene set and to evaluate whether the derived risk score and clinical characteristics serve as independent prognostic factors for overall survival (OS). The prognostic risk score was calculated as ∑(Coef_i × Exp_i), where Coefᵢ denotes the risk coefficient and Exp_i represents the expression level of each gene. Kaplan-Meier survival analysis, executed with the survival and survminer packages in R, was used to compare survival differences between high-risk and low-risk groups, as well as across different gender and age subgroups, with corresponding validation performed on the GEO dataset. The predictive accuracy of the risk score model was assessed by plotting ROC curves using the pROC package.

Principal component analysis (PCA) of prognostic factors

PCA was carried out using the FactoMineR and factoextra packages to investigate the distinctions between high-risk and low-risk groups within both the training and testing cohorts. Visualizations illustrating the clustering based on risk features were also generated.

Construction of nomograms incorporating risk score features

To enhance the reliability of the findings and account for potential confounding variables, univariate and multivariate Cox regression analyses were applied to evaluate the combined effect of the risk score and clinical variables (including gender, age, and pathological stage). The rms package was then used to construct nomograms predicting 1-, 3-, and 5-year survival probabilities for patients. Calibration curves were plotted to assess the concordance between the model-predicted survival rates and the actual observed survival outcomes.

Functional enrichment analysis

To elucidate the potential molecular mechanisms and key biological characteristics underlying the ERG prognostic model, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the clusterProfiler and org.Hs.eg.db packages. Additionally, Gene Set Variation Analysis (GSVA) was conducted to compare biological pathway activities between the high-risk and low-risk groups. The results of these enrichment analyses were visualized with the enrichplot, ggplot2, and pheatmap packages. The hallmark gene sets used for GSVA were sourced from the Molecular Signatures Database (MSigDB).

Analysis of the tumor immune microenvironment

To evaluate disparities in tumor immune cell infiltration between the risk groups, the infiltration scores of 22 immune cell types were estimated using the CIBERSORT algorithm. The results were presented using stacked bar charts and comparative box plots. Furthermore, the expression profiles of genes associated with tumor immune checkpoints were visualized in R using the pheatmap and ggplot2 packages. This analysis aimed to provide deeper insights into the prognostic and immunotherapeutic implications of the tumor immune microenvironment. To characterize the composition of the tumor microenvironment (TME), we employed the ESTIMATE algorithm to compute StromalScore, ImmuneScore, and ESTIMATEScore.

Correlation between prognostic features and TME

We utilized the “survival”, “survminer”, “ggplot2”, and “ggpubr” R packages to compare TME scores between high-risk and low-risk patient groups. To further assess the predictive capacity of the identified genes for immunotherapy response, we evaluated patient prognosis by integrating TME scores with the established risk stratification.

Spatial localization of model genes from bulk RNA-seq via single-cell analysis

Single-cell RNA sequencing (scRNA-seq) data were processed and analyzed using the Seurat package. Initial quality control involved filtering cells based on the following criteria: (I) genes expressed in fewer than 3 cells (to mitigate doublets); (II) total unique molecular identifier (UMI) counts (nCount) below 200 or above 5,000; (III) number of detected genes (nFeature) below 200; (IV) mitochondrial gene content exceeding 20%; and (V) gene count-to-UMI ratio (complexity) below 0.8. Following filtering, gene expression data were normalized and scaled using the NormalizeData and ScaleData functions. To account for cell cycle effects, data transformation was performed using the SCTransform function. PCA was conducted for dimensionality reduction. To correct for batch effects, we applied the Harmony algorithm, computing the top 30 principal components. Cell clustering was then performed using the FindNeighbors and FindClusters functions (with resolution set to 0.6), and results were visualized via Uniform Manifold Approximation and Projection (UMAP) using RunUMAP. Cell type annotation was performed by downloading reference datasets and employing the SingleR package. Model genes identified from bulk RNA-seq were mapped onto the single-cell clusters and visualized using the FeaturePlot function. Expression levels of prognosis-associated differentially expressed genes were quantified across different cell types. To evaluate EMT pathway activity, we applied the AUCell algorithm. This method ranks gene expression within each cell, constructs a ROC curve for a given gene set, and calculates the area under the curve (AUC) as an activity score. The AUCell_buildRankings function was used to establish expression rankings, and AUCell_calcAUC was used to compute scores. Cells were subsequently dichotomized into high and low AUC score groups based on the median value, and their distribution was visualized on UMAP plots.

Cell-cell communication analysis using CellChat

To investigate intercellular signaling networks within the tissue, we performed cell communication analysis on the annotated scRNA-seq data using CellChat (version 1.6.0). The workflow was as follows: first, using the cell type annotations derived from Seurat analysis [including endothelial cells, epithelial cells, monocytes, T cells, B cells, natural killer (NK) cells, and tissue stem cells], we extracted normalized expression matrices and imported them into CellChat to create a communication object. Leveraging CellChat’s curated ligand-receptor interaction database, we inferred potential signaling events between cell populations. The aggregate interaction strength of all ligand-receptor pairs was summarized, and a global communication network was plotted to illustrate the quantity of interactions among different cell types. We further extracted ligands and receptors associated with specific pathways to analyze their communication patterns, identifying sender, receiver, mediator, and influencer cells. By examining outgoing signaling patterns, we elucidated the distinct roles various cell types play in signal propagation. All network visualizations were generated using CellChat’s built-in plotting functions.

Statistical analysis

All statistical analyses were performed using R software (version 4.3.3). Differential expression analysis was conducted using the DESeq2 package with thresholds of |log2FC| >1.5 and P<0.05. Survival analyses were performed using Kaplan-Meier curves with log-rank tests, and Cox proportional hazards regression models were used for univariate and multivariate analyses. LASSO-Cox regression was implemented using the glmnet package with 10-fold cross-validation to select optimal penalty parameters. Receiver operating characteristic (ROC) curves were generated using the pROC package, and the AUC was calculated to evaluate predictive performance. A two-sided P value <0.05 was considered statistically significant.

Results

Characteristic and outcome of the participants

Due to the use of public databases, the available patient information is limited. The patient information obtained in this study is as follows:

The demographic and clinical characteristics of patients in both cohorts are summarized in Table 1. In the TCGA training cohort (n=128), there were 83 males (64.8%) and 45 females (35.2%), with 59 patients (46.1%) aged <60 years and 69 patients (53.9%) aged 60–87 years. Regarding tumor stage, 44 patients (34.4%) were classified as stage I/II and 84 patients (65.6%) as stage III/IV.

In the GSE41613 validation cohort (n=97), there were 66 males (68.0%) and 31 females (32.0%), with 50 patients (51.5%) aged <60 years and 47 patients (48.5%) aged 60–88 years. For tumor stage, 41 patients (42.3%) were stage I/II and 56 patients (57.7%) were stage III/IV. No patients had missing data for these variables. A summary of cohort characteristics, sample selection, and follow-up duration is provided in Table 2.

Differential expression analysis and identification of EMT-associated genes in normal vs. TSCC tissues

Comparative analysis between normal tongue tissues and TSCC samples identified a total of 1,373 DEGs, defined by thresholds of |log₂FC| >1.5 and P<0.05. The distribution of these DEGs is visualized in a volcano plot (Figure 2A). A heatmap was generated to display the expression patterns of the top 25 upregulated and downregulated DEEGs in tumor tissues (Figure 2B). To further illustrate inter-group differences, the ten genes exhibiting the greatest absolute log₂FC values were selected for comparative visualization, comprising seven upregulated and three downregulated genes (Figure 2C). Subsequently, 4,526 ERGs were retrieved from the GeneCards database. Intersection of these ERGs with the identified DEGs yielded 438 DEEGs (Figure 2D).

Figure 2 DEGs screening in tongue cancer based on TCGA database. (A) Volcano plot of DEGs in tongue cancer. (B) Heatmaps of the top 25 upregulated and the top 25 downregulated DEGs. (C) Differential expression of the top 10 DEGs between normal tissues and tongue cancer tissues. (D) Venn diagram showing the overlap of DEGs. *, P<0.05; ***, P<0.001. DEGs, differentially expressed genes; EMT, epithelial-mesenchymal transition; OTC, oral tongue cancer; TCGA, The Cancer Genome Atlas.

Development of a prognostic risk model

Clinical data for TSCC patients were sourced from the TCGA-head and neck squamous cell carcinoma (HNSC) cohort, and after excluding individuals with missing survival information, a total of 128 patients were included in the study, among whom 50 outcome events (deaths) occurred during follow-up. From the 438 DEEGs identified, univariate Cox regression analysis identified 27 genes significantly associated with TSCC prognosis (selecting factors with P value <0.05), as shown in Table 3. To mitigate overfitting during feature selection, LASSO regression analysis was applied, which refined the list to 16 significant genes, with the trajectory of regression coefficients illustrated (Figure 3A). Subsequent multivariate Cox regression analysis on these prognostic genes confirmed the selection of these 16 key genes (Figure 3B), which are potentially implicated in patient outcomes. These analytical steps aid in elucidating prognosis-relevant genes, offering valuable insights for both clinical and fundamental research. A risk score for each patient was calculated based on the expression levels of these genes and their corresponding regression coefficients from the multivariate Cox model, as defined by the following formula: risk score = (−0.113 × MMP13 expression) + (0.458 × HOXA1 expression) + (−0.438 × GRIA3 expression) + (−0.280 × F2RL2 expression) + (0.138 × LINGO2 expression) + (0.359 × ANO1 expression) + (−0.492 × ANO1 expression) + (0.073 × DKK1 expression) + (0.036 × PLK1 expression) + (0.254 × IRX4 expression) + (0.537 × CADM2 expression) + (0.531 × MT3 expression) + (0.165 × KLK3 expression) + (−0.050 × CLDN17 expression) + (−0.033 × MUC4 expression) + (1.058 × MPO expression).

Figure 3 Screening of prognosis-related genes. (A,B) Prognostic genes associated with DEEGs identified by LASSO regression analysis. (C) Waterfall plot of tumor somatic mutations. (D) Chromosomal locations of 16 DEEGs. (E) Heatmap showing correlations among the expression of 16 DEEGs. *, P<0.05; ***, P<0.001. DEEGs, EMT-related differentially expressed genes; EMT, epithelial-mesenchymal transition; LASSO, least absolute shrinkage and selection operator; TMB, tumor mutation burden.

We further examined the somatic mutation landscape of these 16 DEEGs. Among 510 samples analyzed, 66 (12.94%) harbored mutations in these genes, predominantly missense mutations. Analysis of somatic copy number variations (CNVs) revealed the frequency of CNV events across these DEEGs (Figure 3C), and their chromosomal locations are mapped in Figure 3D. Correlation analysis among the 16 DEEGs indicated a generally high degree of co-expression consistency among these ERGs, with the strongest correlations observed between MPO and CADM2, as well as between MPO and F2RL2 (Figure 3E).

Predictive efficacy of the risk signature

To evaluate the model’s capacity to stratify patients and validate its predictive utility, all patients from the TCGA cohort (training set, n=128) and the GSE41613 cohort (validation set, which includes 97 patients, all of whom were included in the study) were dichotomized into high-risk and low-risk groups based on the median risk score. The results demonstrated a positive association between higher risk scores and an increased number of patient deaths. Furthermore, the expression profiles of the 16 signature genes showed significant disparities between the high-risk and low-risk groups in both datasets (Figure 4A,4B). These findings suggest that the risk score model effectively categorizes TSCC patients into distinct prognostic subgroups. PCA can further elucidate the model’s practicality, with results indicating that the EMT related genes involved in model construction can effectively distinguish between high-risk and low risk patient groups, further confirming the model’s accuracy (Figure 4C,4D).

Figure 4 Assessment and validation of the prognostic value of risk scores in the training and test sets. (A,B) Comparison of risk scores between the high- and low-risk groups. Distribution of survival and death in patients from the training set (tongue cancer) and test set (oral squamous cell carcinoma). Expression heatmaps of 16 DEGs in the training and test sets, respectively. (C,D) PCA of high- and low-risk groups. DEGs, differentially expressed genes; OS, overall survival; PCA, principal component analysis.

Independent prognostic significance of risk features

Kaplan-Meier survival analysis was employed to evaluate OS disparities between the high-risk and low-risk groups. The findings revealed that patients in the low-risk group exhibited significantly superior OS compared to those in the high-risk group, a result consistent across both the TCGA training cohort and the GSE41613 validation cohort (P<0.05; Figure 5A,5B). To further investigate the prognostic value of the risk model across diverse patient populations, we conducted subgroup analyses stratified by age and gender. These analyses consistently demonstrated worse survival outcomes for high-risk patients within each subgroup (Figure 5C-5F). Consequently, the established risk model possesses robust prognostic utility for various demographic subgroups of TSCC patients.

Figure 5 Subgroup survival analysis of a prognostic model for DEEGs. (A,B) Kaplan-Meier analysis of risk models on the training and test sets. (C,D) Kaplan-Meier analysis of risk models on age <60 years and age ≥60 years groups. (E,F) Kaplan-Meier analysis of risk models on male and female groups. DEEGs, EMT-related differentially expressed genes; EMT, epithelial-mesenchymal transition; HR, hazard ratio.

Prognostic assessment within the TCGA cohort

To ascertain whether the prognostic features function as independent predictors of survival, univariate and multivariate Cox proportional hazards regression analyses were performed. Univariate analysis identified both clinical stage and the risk score as significant prognostic factors. However, in the subsequent multivariate analysis, only clinical stage retained its independent prognostic significance (Figure 6A,B). To enhance prognostic precision, we integrated the risk score with clinical stage to develop a prognostic nomogram using the TCGA dataset. This nomogram facilitates the prediction of 1-, 3-, and 5-year OS probabilities (Figure 6C). The discriminatory power of the risk features was evaluated using time-dependent ROC curves. The AUC for OS prediction was 0.844, which substantially outperformed the predictive accuracy of age (AUC =0.491), gender (AUC =0.515), and clinical stage alone (AUC =0.630) (Figure 6D). These results underscore the strong predictive capability of the risk model. Furthermore, the model maintained high predictive accuracy for 1-, 3-, and 5-year survival rates, with all corresponding AUC values exceeding 0.75: the 1-year AUC was 0.894 [95% confidence interval (CI): 0.828–0.960], the 3-year AUC was 0.805 (95% CI: 0.702–0.908), and the 5-year AUC was 0.795 (95% CI: 0.643–0.947) (Figure 6E). Calibration curves were plotted to assess the agreement between predicted and observed survival probabilities. The close alignment of these curves with the ideal 45-degree diagonal line confirmed the nomogram’s satisfactory calibration for predicting survival at 1, 3, and 5 years (Figure 6F).

Figure 6 Nomogram construction and independent validation based on prognostic features of ERGs. (A) Univariate Cox regression; (B) Multivariate Cox regression. (C) Nomogram predicting 1-, 3-, and 5-year overall survival in TCGA tongue cancer patients using pathologic stage and risk score. (D) ROC curves comparing risk score and clinical characteristics. (E) ROC curves for 1-, 3-, and 5-year survival predictions. (F) Calibration curves. AUC, area under the curve; CI, confidence interval; EMT, epithelial-mesenchymal transition; ERGs, EMT-related genes; FPR, false positive rate; HR, hazard ratio; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; TPR, true positive rate.

Functional enrichment and gene set enrichment analysis (GSEA)

To elucidate the underlying biological mechanisms distinguishing high-risk and low-risk patients, we performed functional enrichment analyses on differentially expressed genes. GO analysis revealed significant enrichment in biological processes (BPs) related to muscle system processes, muscle contraction, and the regulation of heart contraction. Cellular component (CC) analysis highlighted enrichment in contractile fibers, myofibrils, sarcomeres, and fibrillar collagen trimers. Molecular function (MF) analysis identified enrichment for functions such as extracellular matrix (ECM) structural constituents conferring tensile strength, structural constituents of muscle, platelet-derived growth factor binding, and ligand-gated channel activity (Figure 7A). Subsequently, GSVA was conducted using the R package to perform KEGG pathway analysis. The results, visualized via heatmap, indicated that the high-risk group was characterized by the enrichment of pathways including the Fanconi anemia pathway, translesion synthesis (TLS) and double-strand break (DSB) formation for lesion bypass, and DNA replication licensing. Conversely, the low-risk group showed enrichment in pathways such as the Wnt signaling pathway (Figure 7B).

Figure 7 Enrichment analysis. (A) GO enrichment analysis between the two risk groups. (B) GSVA analysis revealing signaling pathways between high-risk and low-risk groups. (C,D) GSEA analysis between the two different risk groups. BP, biological process; CC, cellular component; GO, Gene Ontology; GSEA, gene set enrichment analysis; GSVA, Gene Set Variation Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Signaling pathways involving Wnt inhibition, the pathogen Yersinia yopH targeting the T cell receptor/nuclear factor of activated T-cells (TCR/NFAT) pathway, and the IL-10 family linked to the JAK/STAT pathway were identified. GSEA revealed that differentially expressed genes in the high-risk group were significantly enriched in pathways associated with dilated cardiomyopathy, drug metabolism via cytochrome P450, oxidative phosphorylation, and retinol metabolism. Conversely, in the low-risk group, differentially expressed genes were significantly enriched in pathways related to cytokine-cytokine receptor interactions, ECM receptor interactions, focal adhesion, and the T cell receptor signaling pathway (Figure 7C,7D).

Analysis of the immune microenvironment based on prognostic risk features

The heterogeneity of the tumor immune microenvironment influences patient prognosis and the efficacy of immunotherapy. In the TME of TSCC, genes such as HLA, MICB, ICAM-1, ICOSLG, and HMGB1 were found to be activated, suggesting that TSCC may be responsive to immune checkpoint inhibitors (Figure 8A). To further investigate the relationship between risk scores and the TME, we employed the CIBERSORT algorithm to estimate the abundance of immune cell infiltration in TSCC patients. Our findings indicated that, compared to the low-risk group, the proportions of activated dendritic cells and activated mast cells showed an increasing trend (Figure 8B,8C). Previous research has established a positive correlation between tumor mutation burden (TMB) and tumor stage, grade, and immune cell infiltration. Using the median TMB value as a cutoff, we stratified TSCC patients into “high TMB” and “low TMB” groups and performed survival analysis. Kaplan-Meier analysis demonstrated that OS was significantly better in the low TMB group compared to the high TMB group. A combined survival analysis incorporating both TMB and tumor risk scores further revealed that patients with low-risk scores and low tumor burden had the most favorable prognosis, whereas those with high-risk scores and high tumor burden had the poorest prognosis (Figure 8D,8E).

Figure 8 Tumor immune microenvironment landscape for ERG models. (A) Expression of immune checkpoint-related genes. (B,C) Different types of immune cells infiltrating tumor tissue in different risk groups. (D) Survival analysis curves of the high-TMB and low-TMB groups. (E) Effect of TMB combined with different risk groups on the probability of survival. *, P<0.05; **, P<0.01. Cor, correlation coefficient; EMT, epithelial-mesenchymal transition; ERG, EMT-related gene; NK, natural killer; ns, no significance; TMB, tumor mutation burden.

Single-cell transcriptome analysis reveals cellular heterogeneity and EMT state characteristics

This study conducted an in-depth analysis of samples using single-cell RNA sequencing data, systematically elucidating the composition of cell types, gene expression profiles, and their association with EMT status. The quality control assessment (Figure 9A-9C) demonstrated high consistency across samples in terms of the number of RNA molecules, detected genes, gene-to-UMI count ratios, and the proportion of mitochondrial genes, confirming data reliability for subsequent analyses. Cell type identification results (Figure 9D) indicated the presence of various cell populations, including tissue stem cells, endothelial cells, NK cells, T cells, epithelial cells, and monocytes. In the gene expression visualization (Figure 9E), MMP13 and ANO1 displayed distinct regional expression patterns on the UMAP plot, suggesting their potential roles in regulating specific cellular subpopulations. Heatmap analysis (Figure 9F) further illustrated the expression patterns of marker genes across different cell types, highlighting transcriptional heterogeneity. Assessment of EMT activity (Figure 9G) showed that epithelial cells and tissue stem cells exhibited higher EMT_AUC scores, indicating a greater potential for undergoing EMT. Finally, cells were categorized into high-EMT and low-EMT groups based on their EMT scores (Figure 9H), providing a clear framework for subsequent investigations into the role of EMT in cellular state transitions and functional regulation. These findings elucidate the relationship between cellular heterogeneity and EMT dynamics in TSCC.

Figure 9 Single-cell RNA sequencing analysis reveals cellular heterogeneity and EMT-related gene expression in tumor samples. (A-C) QC metrics for single-cell data: (A) pre-QC (UMI, genes, complexity, mito%); (B) post-QC; (C) post-QC correlation (UMI vs. genes, R²=0.99). QC thresholds: complexity >0.8, mito% <20%. (D) UMAP plot showing cell types identified in the sample. (E) Expression of MMP13, ANO1 and PLK1 across the UMAP space. (F) Average expression of selected genes in each cell type. Dot size = percentage of expressing cells; color = expression level. (G) EMT activity score (EMT_AUC) across cell types. (H) UMAP plot with cells grouped into high or low EMT activity. AUC, area under the curve; CMP, common myeloid progenitor; EMT, epithelial-mesenchymal transition; NK, natural killer; QC, quality control; UMAP, Uniform Manifold Approximation and Projection; UMI, unique molecular identifier.

Heterogeneity, key gene expression, and EMT status at the single-cell level: providing critical insights into cellular state dynamics in physiological and pathological contexts

This research systematically investigated the intricate signaling communication networks among diverse cell populations within the tissue microenvironment using advanced cell communication analysis techniques. Initially, the study identified the key cell types participating in essential communications (Figure 10A), encompassing endothelial cells, epithelial cells, monocytes, T cells, B cells, NK cells, and tissue stem cells, thereby establishing a foundation for subsequent network analysis. By quantifying the frequency of interactions between these cell types (Figure 10B), the analysis revealed significant and recurrent communication between epithelial cells and immune cells, as well as between tissue stem cells and adjacent stromal cells. This observation suggests these interactions may constitute a central mechanism in maintaining microenvironmental homeostasis or in driving pathological progression. Focusing further on critical pathways involved in ECM remodeling, the study delineated a detailed signaling pathway network (Figure 10C). This network clarified the specific roles of signaling senders, receivers, mediators, and regulators within this pathway, highlighting the pivotal function of collagen signaling in mediating interactions between cells and the stroma and in regulating cellular behavior. Finally, through a quantitative assessment of outgoing signaling patterns (Figure 10D), the relative signaling output strength of each cell type was measured. The results indicated that tissue stem cells and specific immune cell subsets demonstrate particularly active signaling output, potentially acting as primary drivers of the overall microenvironmental communication network. In conclusion, this study elucidates, from a systems-level perspective, a cell communication architecture predominantly orchestrated by collagen signaling pathways. It clarifies the distinct contributions of various cell populations within these signaling networks, offering valuable insights into network dynamics and a molecular foundation for understanding cooperative cellular mechanisms in tissue development, homeostasis, and the onset and progression of disease.

Figure 10 Intercellular communication analysis reveals collagen signaling network and cellular interaction patterns in the TSCC microenvironment. (A) Intercellular communication analysis reveals collagen signaling network and cellular interaction patterns. (B) Histogram displaying the total number of inferred interactions between different cell types. (C) Schematic network diagram of the “COLLAGEN signaling pathway”. Nodes are labeled according to their functional roles in communication: Sender, Receiver, Mediator, and Influencer. (D) Heatmap showing the outgoing signaling patterns of each cell type across different signaling pathways. Rows represent cell types; columns represent pathways. Color indicates the strength of outgoing signals (scaled). Bar plots show the relative strength of each cell type in the sample. (E) Cell-cell communication probability of ligand-receptor pairs among tumor microenvironment cell subsets. CMP, common myeloid progenitor; NK, natural killer.

Discussion

TSCC is a prevalent malignant neoplasm of the oral cavity, with its incidence increasing in numerous regions globally (14,15). As a highly aggressive subtype of head and neck squamous cell carcinoma, TSCC presents a substantial threat to patient survival and quality of life. The prognosis for TSCC patients remains unfavorable, largely attributable to high rates of metastasis and recurrence, compounded by challenges related to treatment resistance and the significant adverse effects associated with conventional therapies such as surgery, chemotherapy, and radiotherapy (16,17). The identification of novel prognostic biomarkers and therapeutic targets has become increasingly critical, as current clinical strategies frequently fail to address the complex biological behaviors of this malignancy, leading to suboptimal patient outcomes (8,18).

The EMT is a dynamic and reversible BP strongly associated with cancer metastasis, therapeutic resistance, and poor clinical prognosis across various cancer types (19,20). Recent progress in multi-omics technologies has highlighted the complexity of EMT, demonstrating that it is not a simple binary switch but rather encompasses a spectrum of intermediate cellular states, each with distinct molecular characteristics and clinical relevance (21,22). In this investigation, we employed integrated bioinformatics approaches and single-cell resolution data to develop and validate a novel EMT-related prognostic signature for TSCC. Our 16-gene signature effectively characterizes the aggressive, mesenchymal-like phenotype associated with poor prognosis.

Our analysis began by identifying DEGs from TCGA data and cross-referencing them with a selected EMT gene set, a standard method for extracting functionally significant signatures. Subsequent univariate Cox and LASSO-Cox regression analyses refined the selection to 16 core genes with independent prognostic significance. This data-driven method for signature development, similar to recent pan-cancer and cancer-specific studies, reduces bias and increases biological relevance. The derived risk score demonstrated strong independent prognostic value, surpassing conventional clinicopathological factors such as disease stage and age in time-dependent ROC analyses (Figure 6D,6E). This supports the emerging view that molecular signatures capturing essential BPs, like EMT, often yield better predictive accuracy for patient outcomes than traditional parameters alone (23,24).

The performance of our 16-gene EMT-related signature was rigorously evaluated through external validation using the independent GSE41613 cohort. The development of this signature combined two datasets: the TCGA-TSCC cohort (n=128) for model derivation and the GSE41613 cohort (n=97) for external validation. In the TCGA training cohort, the risk score successfully stratified patients into high- and low-risk groups with significantly divergent OS (P<0.05, Figure 5A), yielding an AUC of 0.844 for OS prediction (Figure 6D), with 1-, 3-, and 5-year AUCs of 0.894, 0.805, and 0.795, respectively. When applied to the GSE41613 validation cohort, the same risk score formula demonstrated consistent prognostic stratification capability, with Kaplan-Meier analysis confirming significantly worse outcomes for high-risk patients (P<0.05, Figure 5B). The predictive accuracy remained robust in the validation set, with time-dependent AUC values of 0.894, 0.805, and 0.795 for 1-, 3-, and 5-year survival, respectively (Figure 6E), closely mirroring the performance metrics obtained from the training cohort.

Several studies have previously developed prognostic signatures for OSCC using various gene sets. Ai and colleagues developed a 9-gene EMT-related signature for OSCC using TCGA data, reporting successful external validation in the GSE41613 dataset with well-performing ROC curves for 1-, 3-, and 5-year survival. Their study also confirmed that AREG, COL5A3, DKK1, GAS1, GPX7 and PLOD2 were distinctly upregulated and SFRP1 downregulated in OSCC compared to normal tissues (25). Fang and colleagues constructed an invasion-related 6-gene prognostic model (HMGN2, MYL12B, ACTB, PPP1CA, PSMB9, IFITM3) for TSCC using TCGA and GSE41116 datasets, demonstrating that patients in the low-risk group had longer disease-free survival than those in the high-risk group (26). Chen and colleagues recently developed a three-gene signature (CXCL12, PLAU, PXDN) for OSCC using multi-algorithm approaches including LASSO and Random Forest, reporting 1/3/5-year AUCs of 0.767/0.625/0.714 in the training cohort and validating across three independent external cohorts with consistent prognostic performance (27). Compared with these existing models, our 16-gene signature demonstrates comparable or superior predictive performance, with consistently high AUC values in both development and validation cohorts. The systematic review by Giles and Rothwell (28), although focused on stroke risk prediction, highlights the importance of external validation in demonstrating model generalizability—they found a pooled AUC of 0.72 (95% CI: 0.63–0.82) for all studies meeting their search criteria, and an AUC of 0.69 (95% CI: 0.64–0.74) after excluding the original derivation studies. Our signature’s performance in both development and validation cohorts exceeds these pooled estimates, suggesting robust predictive capability.

Notably, the consistent performance across two independent datasets—TCGA (primarily North American population) and GSE41613 (mixed population)—supports the generalizability of our signature across different demographic and clinical settings. The slightly higher 1-year AUC in the validation cohort (0.894 vs. 0.844 overall) may reflect differences in case mix, follow-up duration, or population characteristics between the two cohorts. Nevertheless, the narrow range of AUC values across different time points and datasets indicates stable and reproducible predictive accuracy, reinforcing the clinical utility of our EMT-related signature.

Functional enrichment analysis of DEGs between high- and low-risk groups provided robust validation of the biological foundation of our signature. The high-risk group showed significant enrichment for terms related to muscle contraction, ECM organization, and cell adhesion (Figure 7A). This pattern is indicative of a typical mesenchymal state, characterized by increased cytoskeletal dynamics and active ECM remodeling—processes vital for cell migration, invasion, and the formation of a tumor-promoting niche (29,30). The emphasis on ECM-related pathways is especially notable, as recent research has shown that a remodeled, stiffened ECM not only promotes invasion but also triggers mechanosignaling pathways that further enhance EMT and stemness (31-33). Thus, our signature represents a tumor phenotype with high metastatic potential and inherent resilience.

A key aspect of our study was elucidating the immune landscape linked to the EMT-high risk group. Analysis of immune cell infiltration revealed significant differences between risk groups. Longitudinal survival data indicated notably worse outcomes associated with specific immune subsets in certain contexts (Figure 8), suggesting a dysfunctional or excluded immune response in high-risk patients. This finding aligns with extensive recent literature describing the immunosuppressive characteristics of mesenchymal tumor states (34). EMT can promote the secretion of immunosuppressive cytokines (e.g., TGF-β), upregulate checkpoint ligands (e.g., PD-L1), and recruit regulatory immune cells, collectively creating an immune-cold TME that is resistant to checkpoint blockade immunotherapy. Therefore, our signature may be useful in identifying patients less likely to respond to immunotherapies, enabling more personalized treatment approaches.

To overcome the limitations of bulk sequencing and explore cellular heterogeneity, we conducted scRNA-seq. This allowed clear identification of major CCs within the TME, including epithelial, stromal, and immune lineages (Figure 9D). Using gene set scoring methods (e.g., AUCell), we mapped EMT activity at single-cell resolution, confirming that high-EMT programs are restricted to specific subpopulations rather than being uniformly expressed (Figure 9G,9H) (35). This heterogeneity is a crucial aspect of tumor plasticity, where subpopulations with varying EMT states collaborate to drive progression and therapy resistance (36). The distinct expression patterns of signature genes such as MMP13 across clusters (Figure 9E,9F) further identify the cellular players involved in ECM degradation and invasion.

Using CellChat analysis, we modeled the intercellular communication networks rewired in the context of high EMT activity (Figure 10) (37). The predicted changes in signaling intensity and patterns, particularly within collagen and other ECM-related pathways, indicate a self-perpetuating cycle of stromal-epithelial interaction. In this proposed model, EMT-active tumor cells and activated cancer-associated fibroblasts (CAFs) are likely involved in bidirectional signaling through ECM components and secreted factors. This interaction results in progressive fibrosis, matrix stiffening, and the persistent activation of pro-survival and pro-invasive signaling pathways in tumor cells (38,39). This systems-level perspective highlights that the poor prognosis linked to our signature stems not only from intrinsic alterations in tumor cells but also from the co-evolution of a supportive and interactive TME.

To transform our findings into a clinically useful tool, we constructed a nomogram that integrates the risk score with key clinical variables (Figure 6C). Such combined models are gaining recognition for their potential to enhance individualized patient prognosis prediction and aid clinical decision-making (40,41).

This research has several constraints. We acknowledge a limitation of this study: our findings are derived exclusively from public datasets (TCGA and GEO), which may introduce population and technical biases. Recent studies have confirmed the presence of site-specific bias in TCGA data, where deep learning models can classify data acquisition sites with high accuracy based solely on embedded features, suggesting that reliance on biased features may lead to over-optimistic performance estimates (42). Additionally, expression-level dependent biases have been identified in TCGA RNA-seq data that persist after conventional normalization and can corrupt gene-gene correlation estimations (43). However, this limitation is not unique to our study but reflects current practice in the field. Multiple recent prognostic signature studies have successfully utilized the same TCGA + GEO framework for discovery and validation. For instance, a 2025 study developed an 8-gene EMT signature for lung adenocarcinoma using TCGA as training and GSE30219 for validation, acknowledging that “our study used public databases” and that further validation is needed (44). Similarly, a 15-gene EMT signature for breast cancer was validated across 4 independent datasets using the TCGA + GEO framework (45). A recent lung adenocarcinoma study also employed single-cell RNA-seq with TCGA and GEO bulk data to construct prognostic models (46). To mitigate bias, we performed external validation using an independent GEO cohort (GSE41613) with cross-platform consistency (RNA-seq vs. microarray). The risk score successfully stratified patients in the validation cohort (Figure 5B), and we have transparently reported demographic characteristics of both cohorts (Table 1). Prospective multi-center validation remains essential.

An important finding requires explicit discussion: the risk score did not retain independent prognostic significance when clinical stage was included in the model (Figure 6B). This indicates that the prognostic information captured by the signature is partially correlated with or mediated through conventional clinicopathological factors. This phenomenon has been observed in other studies as well. For example, a recent lung adenocarcinoma study that integrated single-cell and bulk transcriptomic data to construct an EMT-associated CAF prognostic model reported that the risk score effectively stratified patients but noted that combining risk scores with clinical variables achieved the highest predictive accuracy (46). Similarly, a 2026 study in Cancer Gene Therapy that reconstructed the EMT continuum in lung adenocarcinoma demonstrated that while ERG signatures are strongly associated with prognosis, their optimal clinical utility is achieved when integrated with clinical parameters (22). Therefore, we propose that our 16‑gene EMT signature is best used in combination with TNM staging—as in our nomogram (Figure 6C)—rather than as a replacement for clinical assessment. This “additive value” approach is consistent with recently published EMT‑based prognostic models and represents a clinically pragmatic strategy.

Additionally, the proposed biological mechanisms remain speculative due to the absence of in vitro and in vivo experimental validation. Recent OSCC studies provide frameworks for such validation. Mozalbat et al. analyzed SMAD4 and EMT markers in 23 OSCC samples and validated findings in an OSCC cell model with SMAD4 mutation, demonstrating correlation between EMT markers and pathological staging (47). Zhang et al. showed that fisetin inhibits periodontal pathogen-induced EMT in OSCC cells via the Wnt/β-catenin pathway, using qRT-PCR and Western blot validation (48). Future studies should prioritize functional characterization of our 16-gene signature using similar approaches (CRISPR-based knockdown in CAL27/SCC-9 cells, followed by xenograft validation).

Furthermore, although a 16-gene signature may appear complex for routine clinical application, it is feasible for translation into a targeted PCR-based panel. A recent study successfully developed and validated a qPCR-based 10-gene mRNA test for rapid (1-hour) prognostic stratification of nasopharyngeal carcinoma, demonstrating that multigene signatures (10–16 genes) can be clinically deployed using routine qPCR platforms (49). Compared with previously reported TSCC prognostic signatures (3-gene, 6-gene, 9-gene models), our 16-gene model provides competitive and stable AUC values (0.894 for 1-year survival in training cohort), supporting its potential clinical utility.

Although the signature was developed and validated using TCGA data, prospective validation in independent, multi-center cohorts is crucial to confirm its broader applicability. The mechanistic functions of the specific 16 genes within the EMT network of TSCC necessitate functional validation via in vitro and in vivo experiments. Moreover, while the scRNA-seq data are informative, they represent a static snapshot; longitudinal or spatial single-cell analyses could more effectively capture the dynamics of EMT state transitions and cell-cell interactions during tumor progression (50).

Conclusions

In summary, we have developed and extensively characterized a 16-gene EMT-related signature that reliably predicts poor prognosis in TSCC. By integrating bulk transcriptomic, survival, and single-cell data, we show that this signature identifies a tumor subtype characterized by a mesenchymal phenotype, an immunosuppressive and fibrotic TME, and altered cellular communication networks. This signature not only represents a promising prognostic biomarker but also reveals potential therapeutic targets. Future research should concentrate on targeting the identified crosstalk pathways (e.g., collagen signaling) or reversing the EMT state itself, potentially in combination with immunotherapy, to disrupt the supportive ecosystem of high-risk tumors and improve patient outcomes (30).

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0435/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0435/prf

Funding: The work was supported by Medical Science Research Project of Hebei (No. 20230135).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0435/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Jeong S, Choi HI, Yang KI, Kim JS, Ryu JW, Park HJ. Artificial Intelligence in the Diagnosis of Tongue Cancer: A Systematic Review with Meta-Analysis. Biomedicines 2025;13:1849. [Crossref] [PubMed]
2. Gonzalez M, Riera March A. Tongue Cancer. In: StatPearls [Internet]. Treasure Island (FL): StatPearls Publishing; 2026 Jan-. Available online: https://www.ncbi.nlm.nih.gov/books/NBK562324/
3. Qing B, Li X, He X, et al. Characterizing Stage-Specific Cellular Dynamics and Microenvironmental Remodeling in Lung Adenocarcinoma by Single-Cell RNA Sequencing. Adv Sci (Weinh) 2026;13:e10847. [Crossref] [PubMed]
4. Li C, Enciso-Martinez A, Zhu L, et al. Surface-Associated Proteins on Extracellular Vesicles Remodel the Tumor Microenvironment by Potentiating TGF-β Signaling in a Contact-Dependent Manner. Adv Sci (Weinh) 2026;13:e13286. [Crossref] [PubMed]
5. Shen S, Zhou H, Xiao Z, et al. PRMT1 in human neoplasm: cancer biology and potential therapeutic target. Cell Commun Signal 2024;22:102. [Crossref] [PubMed]
6. Sun ED, Ma R, Zou J. SPRITE: improving spatial gene expression imputation with gene and cell networks. Bioinformatics 2024;40:i521-8. [Crossref] [PubMed]
7. Haupt L. Authenticity and Clinical Decision-Making. Hastings Cent Rep 2022;52:2. [Crossref] [PubMed]
8. Zang Z, Xu Y, Lu L, et al. UDRN: Unified Dimensional Reduction Neural Network for feature selection and feature projection. Neural Netw 2023;161:626-37. [Crossref] [PubMed]
9. Lin W, Huang W, Mei S, et al. Transdifferentiation of tongue muscle cells into cancer-associated fibroblasts in response to tongue squamous cell carcinoma. Nat Commun 2025;16:6753. [Crossref] [PubMed]
10. Liu XF, Wang B, Zhou M, et al. Integrative analysis of EMT-driving genes identifies a prognostic signature and GJB2 as a potential biomarker in glioblastoma. Front Cell Dev Biol 2025;13:1754988. [Crossref] [PubMed]
11. Yang H, Zhu H, Ahn M, et al. Weighted functional linear Cox regression model. Stat Methods Med Res 2021;30:1917-31. [Crossref] [PubMed]
12. Murphy S. Principles of Tumor Biology. Vet Clin North Am Equine Pract 2024;40:341-50. [Crossref] [PubMed]
13. Xie Y, Shi H, Han B. Bioinformatic analysis of underlying mechanisms of Kawasaki disease via Weighted Gene Correlation Network Analysis (WGCNA) and the Least Absolute Shrinkage and Selection Operator method (LASSO) regression model. BMC Pediatr 2023;23:90. [Crossref] [PubMed]
14. Almangush A, Salo T, Haglund C, et al. Validation of a New Histopathologic Risk Model in Early Oral Tongue Cancer: A Combination of a Modified Worst Pattern of Invasion and a New Tumor Budding Score. Am J Surg Pathol 2025;49:1036-41. [Crossref] [PubMed]
15. Burus T, Damgacioglu H, Huang B, et al. Trends in Oral Tongue Cancer Incidence in the US. JAMA Otolaryngol Head Neck Surg 2024;150:436-43. [Crossref] [PubMed]
16. Engle JA, Dibb JT, Jakob JA. High Tumor Mutational Burden in Hepatocellular Carcinoma. Cureus 2024;16:e68132. [Crossref] [PubMed]
17. Schaafsma E, Zhao Y, Wang Y, et al. Whole transcriptome signature for prognostic prediction (WTSPP): application of whole transcriptome signature for prognostic prediction in cancer. Lab Invest 2020;100:1356-66. [Crossref] [PubMed]
18. Zhong L, Yang F, Sun S, et al. Predicting lung cancer survival prognosis based on the conditional survival bayesian network. BMC Med Res Methodol 2024;24:16. [Crossref] [PubMed]
19. Yang J, Antin P, Berx G, et al. Guidelines and definitions for research on epithelial-mesenchymal transition. Nat Rev Mol Cell Biol 2020;21:341-52. [Crossref] [PubMed]
20. Lambert AW, Weinberg RA. Linking EMT programmes to normal and neoplastic epithelial stem cells. Nat Rev Cancer 2021;21:325-38. [Crossref] [PubMed]
21. Pastushenko I, Brisebarre A, Sifrim A, et al. Identification of the tumour transition states occurring during EMT. Nature 2018;556:463-8. [Crossref] [PubMed]
22. Qu Q, Ma Y, Huang C, et al. Construction of the cancer cell continuum reveals hybrid EMT state driving lung adenocarcinoma aggression. Cancer Gene Ther 2026;33:186-97. [Crossref] [PubMed]
23. Horndalsveen H, Haakensen VD, Madebo T, et al. Blood-based tumor mutational burden as a biomarker in unresectable non-small cell lung cancer treated with chemoradiotherapy and durvalumab. Front Oncol 2025;15:1681420. [Crossref] [PubMed]
24. Wang J, Zhang W, Zhang J, et al. Integrative Multi-Omics Analysis Uncovers Immunological Phenotypes Predictive of Combinatorial Immunotherapy Response in Gastric Cancer. Adv Sci (Weinh) 2026;13:e14482. [Crossref] [PubMed]
25. Ai J, Tan Y, Liu B, et al. Systematic establishment and verification of an epithelial-mesenchymal transition gene signature for predicting prognosis of oral squamous cell carcinoma. Front Genet 2023;14:1113137. [Crossref] [PubMed]
26. Fang W, Chen S, Wan D, et al. Identification and Validation of an Invasion-Related Disease-Free Survival Prognostic Model for Tongue Squamous Cell Carcinoma. Oncology 2025;103:237-52. [Crossref] [PubMed]
27. Chen J, Kim D, Kim JY, et al. Development and multi-cohort validation of a prognostic risk score model for oral squamous cell carcinoma based on a three-gene signature. Cancer Genet 2025;298-299:88-98. [Crossref] [PubMed]
28. Giles MF, Rothwell PM. Systematic review and pooled analysis of published and unpublished validations of the ABCD and ABCD2 transient ischemic attack risk scores. Stroke 2010;41:667-73. [Crossref] [PubMed]
29. Sharma A, Steger RF, Li JM, et al. Sp1 mechanotransduction regulates breast cancer cell invasion in engineered viscoelastic extracellular matrices. Biomaterials 2026;327:123755. [Crossref] [PubMed]
30. Tiskratok W, Kyawsoewin M, Thuephut R, et al. Matrix Stiffness Drives Aggressive Phenotype in Tongue Squamous Cell Carcinoma via Mechanotransduction-Stromal Signalling. Int Dent J 2026;76:109581. [Crossref] [PubMed]
31. Boumahdi S, de Sauvage FJ. The great escape: tumour cell plasticity in resistance to targeted therapy. Nat Rev Drug Discov 2020;19:39-56. [Crossref] [PubMed]
32. De Blander H, Marine JC. A Physical Framework to Control Cancer Cell Heterogeneity and Plasticity. Cancer Discov 2026;16:637-43. [Crossref] [PubMed]
33. Kazemi KS, Miyazawa M, Hanemann JAC, et al. A Pan-Cancer Transcriptomic Signature for Conserved Molecular Programs Underlying Premalignant-Malignant Progression Across Common Carcinomas. Dent J (Basel) 2026;14:228. [Crossref] [PubMed]
34. Dongre A, Rashidian M, Eaton EN, et al. Direct and Indirect Regulators of Epithelial-Mesenchymal Transition-Mediated Immunosuppression in Breast Carcinomas. Cancer Discov 2021;11:1286-305. [Crossref] [PubMed]
35. Andreatta M, Carmona SJ. UCell: Robust and scalable single-cell gene signature scoring. Comput Struct Biotechnol J 2021;19:3796-8. [Crossref] [PubMed]
36. Lüönd F, Sugiyama N, Bill R, et al. Distinct contributions of partial and full EMT to breast cancer malignancy. Dev Cell 2021;56:3203-3221.e11. [Crossref] [PubMed]
37. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
38. Li Q, Yang C, Li J, et al. The type I collagen paradox in PDAC progression: microenvironmental protector turned tumor accomplice. J Transl Med 2025;23:744. [Crossref] [PubMed]
39. Song Z, Ren L, Wang H, et al. Single Cell and Bulk RNA-Seq Profiling of Non-Metastatic Versus Bone-Metastatic Prostate Cancer Identifies the CXCL10-CXCR3 Axis as a Key Determinant of Tumor Microenvironment and Treatment Resistance. Biomedicines 2026;14:943. [Crossref] [PubMed]
40. Zeng C, Xu C, Wei Y, et al. Training and experimental validation a novel anoikis- and epithelial‒mesenchymal transition-related signature for evaluating prognosis and predicting immunotherapy efficacy in gastric cancer. J Cancer 2025;16:1078-100. [Crossref] [PubMed]
41. Li Y, Pan W, Ying X, et al. Integrated single-cell and bulk transcriptomic analyses unveil a necroptosis-related prognostic model and its association with tumor microenvironment remodeling in gastric cancer. Discov Oncol 2026; Epub ahead of print. [Crossref]
42. Kheiri F, Rahnamayan S, Makrehchi M, et al. Investigation on potential bias factors in histopathology datasets. Sci Rep 2025;15:11349. [Crossref] [PubMed]
43. Thron C, Jafari F. Correcting scale distortion in RNA sequencing data. BMC Bioinformatics 2025;26:32. [Crossref] [PubMed]
44. Han P, Guda C, Liu Q. An efficient epithelial-mesenchymal transition-related gene signature for predicting the survival of patients with lung adenocarcinoma. Transl Cancer Res 2025;14:4989-5001. [Crossref] [PubMed]
45. Liang W, Wang ZY, Shao QF, et al. Genes From Epithelial-Mesenchymal Transition Predict Overall Survival Effectively in Breast Cancer: A Novel Risk Model Based on Initial Step of Tumor Metastasis. Breast Cancer (Auckl) 2026;20:11782234261433697. [Crossref] [PubMed]
46. An H, An P. Single-cell and bulk RNA analysis identifies EMT-associated CAF signatures and prognostic model in lung adenocarcinoma. Discov Oncol 2025;16:1106. [Crossref] [PubMed]
47. Mozalbat S, Nashef A, Maalouf N, et al. The Interplay of SMAD4 and EMT in Oral Squamous Cell Carcinoma. Cancers (Basel) 2025;17:1761. [Crossref] [PubMed]
48. Zhang R, Takigawa H, Maruyama H, et al. Fisetin Inhibits Periodontal Pathogen-Induced EMT in Oral Squamous Cell Carcinoma via the Wnt/β-Catenin Pathway. Nutrients 2025;17:3522. [Crossref] [PubMed]
49. Liang Y, Mo Z, Teh MT. A Multigene Signature for Prognostic Stratification of Nasopharyngeal Carcinoma. Cancers (Basel) 2026;18:1197. [Crossref] [PubMed]
50. Zhuang X. Spatially resolved single-cell genomics and transcriptomics by imaging. Nat Methods 2021;18:18-22. [Crossref] [PubMed]