Comprehensive analysis and validation of the prognostic value of glycosylation- and mitophagy-related genes for head and neck squamous cell carcinoma

Introduction

An estimated 3% of cancer diagnoses and 1.5% of cancer deaths in the United States are caused by head and neck squamous cell carcinoma (HNSCC), according to the International Agency for Research on Cancer (1). HNSCC is a malignant tumor associated with low survival rates and poor quality of life. Despite surgery, radiation, and chemotherapy being the principal therapies for HNSCC, clinical results are unsatisfactory due to tumor heterogeneity and therapeutic resistance (2).

New research is shedding insights into the tumor microenvironment as well as the molecular pathways that promote HNSCC progression. Among these, glycosylation and mitophagy have emerged as pivotal processes influencing tumor metabolism, immune evasion, and clinical outcomes (3-7). Glycosylation, a fundamental post-translational modification, regulates cell signaling and immune interactions, with aberrant patterns frequently associated with tumor aggressiveness and poor prognosis in HNSCC (8). The process of mitophagy, which controls cell death and energy homeostasis, is additionally associated to cancer cell proliferation and survival (9), particularly in HNSCC (10).

Despite extensive research, the interaction between glycosylation and mitophagy in HNSCC remains inadequately understood. While both processes have been individually linked to cancer progression, their combined impact on prognosis and treatment response has not been thoroughly investigated. Bridging this information gap may enable the discovery of new prognostic biomarkers and treatment targets.

To explore these interactions, we used bioinformatics analysis on The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets. This approach allowed for a comprehensive assessment of glycosylation- and mitophagy-related genes (G&MRGs) markers in HNSCC. By elucidating these molecular relationships, our findings might contribute to early diagnosis and inform personalized treatment strategies, ultimately improving clinical outcomes and patient survival. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0011/rc).

Methods

Data sources

We used the TCGAbiolinks package (11) (Version 2.30.0)in the R software to download data from TCGA database (https://portal.gdc.cancer.gov/), specifically the HNSCC dataset (TCGA-HNSCC). After excluding samples without clinical data, we obtained the count-formatted sequencing data for 503 HNSCC and 44 control samples, which were normalized to the “fragments per kilobase of transcript per million mapped reads” (FPKM) format. Clinical data were retrieved from the UCSC Xena database (12). The relevant details are presented in (Table 1). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All data used in this study were obtained from publicly available databases, and no additional informed consent was required.

We used the GEOquery package (13,14) (Version 2.70.0)in R to download the GSE107591 (15,16) and GSE6631 (17) datasets, which included data on head and neck tissue samples obtained from Homo sapiens. GSE107591 contained 24 HNSCC and 23 control samples, whereas GSE6631 contained 22 HNSCC and 22 control samples. Comprehensive details about the samples are shown in (Table 2). However, the GEO datasets (GSE107591 and GSE6631) contain only limited baseline information and lack detailed clinicopathological annotations. Because the GEO datasets used in this study did not contain complete follow-up information required for survival analysis, external survival validation could not be conducted. To ensure the comparability of gene expression levels across different GEO datasets, the downloaded expression matrices were uniformly processed and integrated into a combined dataset. Batch effects resulting from differences in experimental platforms and processing conditions were corrected using the ComBat() function from the sva (Version 3.50.0) (18) R package. This step was performed to minimize potential systematic biases across datasets. The integrated data were subsequently normalized using the limma (Version 3.58.1) (19) R package. Differential expression analysis was then performed based on linear model fitting followed by empirical Bayes moderation using the eBayes() function.

Gene information was obtained from GeneCards (20), yielding 2,882 glycosylation-related genes (GRGs) and 4,878 mitophagy-related genes (MRGs). A PubMed search using the keywords “Glycosylation” and “Mitophagy” further identified additional 285 GRGs (21) and 29 MRGs (22). After integration, 866 overlapping genes were identified as G&MRGs.

Batch effects were corrected using the ComBat function implemented in the sva R package. To assess the effectiveness of batch effect correction, principal component analysis (PCA) (23) was conducted to visualize the clustering patterns of samples before and after correction.

Glycosylation- and mitophagy-related differentially expressed genes (G&MRDEGs) associated with HNSCC

Based on the sample grouping in the TCGA-HNSCC dataset, the samples were separated into two groups: one for HNSCC and one for control. The R package limma was used to analyze differentially expressed genes (DEGs) in these groups, using the following thresholds: |logFC| >2 and adj. P<0.05. Upregulated genes were those with logFC >2 and adj. P<0.05, whereas downregulated genes were those with logFC <−2 and adj. P<0.05. Correcting P values was done using the Benjamini-Hochberg (BH) approach. Using the ggplot2 (Version 3.4.4) program in R, we created a volcano plot to display the data.

Differential expression analysis for the integrated GEO dataset was performed using the limma R package. The resulting P-values were adjusted using the BH method to control the false discovery rate. Genes with |logFC| >0.25 and adjusted P<0.05 were considered DEGs. Specifically, genes with logFC >0.25 and adjusted P<0.05 were defined as upregulated, whereas genes with logFC <−0.25 and adjusted P<0.05 were defined as downregulated.

We combined G&MRGs from the GEO datasets with DEGs from TCGA-HNSCC and intersected them to discover G&MRDEGs associated with HNSCC. R packages pheatmap (Version 1.0.12) and RCircos (24) (version 1.2.2) were employed to generate heatmaps and circular plots, respectively, in the visualization process.

Somatic mutation (SM) and copy number variation (CNV) analysis

To examine SMs in the HNSCC cohort within the TCGA-HNSCC dataset, we extracted the “masked somatic mutation” data, preprocessed it with VarScan, including formatting mutation information and extracting mutation type annotations. Subsequent analyses and visualization were performed using the maftools (Version 1.42.0) (25) R package. These analyses included classification of mutation types, analysis of single nucleotide variant (SNV) patterns, ranking genes by mutation frequency, and generating oncoplots to visualize the mutational landscape. In the CNV research of the HNSCC cohort, we employed the “masked copy number segment” obtained from TCGA and performed the analysis via GISTIC2.0 (26). This algorithm was used to identify significantly amplified and deleted genomic regions and to estimate the frequency of gene-level copy number alterations across the cohort.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis

GO analysis (27) and the KEGG (28) pathway analysis were performed based on the 23 G&MRDEGs to explore their potential biological functions and associated signaling pathways. The enrichment results were visualized using bubble plots and network diagrams generated with the ggplot2 package and relevant functions from the clusterProfiler (29) (Version 4.10.0) package in R. Raw P values were adjusted for multiple testing using the BH method, and significance levels established at P<0.05 and FDR <0.25. This FDR cutoff was adopted according to the recommendations of the gene set enrichment analysis (GSEA) (30) methodology to balance false-positive control with sensitivity for detecting potential biological signals in exploratory analyses. Pathway visualization was conducted utilizing the R package Pathview (Version 1.42.0) (31).

GSEA

GSEA was performed with the TCGA-HNSCC dataset. Genes were ranked based on their logFC values obtained from the differential expression analysis described in Section 2.2. Enrichment analysis was conducted utilizing the R package clusterProfiler (Version 4.10.0) with the following parameters: seed =2020, permutations =1,000, and gene set size 10-500. The investigation utilized gene sets from the Molecular Signatures Database (MSigDB) (32). The criterion for significance were adj. P<0.05 and FDR <0.25, with the BH method employed for P value adjustment.

Construction of a predictive risk model for HNSCC

A predictive risk model for HNSCC was developed utilizing the TCGA-HNSCC dataset and the R package survival (33) (Version 3.5-7). Univariate Cox regression analysis was performed to screen variables with P<0.10, which were subsequently included in least absolute shrinkage and selection operator (LASSO) regression analysis using the R package glmnet (34) (Version 4.1-8). LASSO regression, which incorporates a penalty term (lambda × absolute slope), minimizes overfitting. The results were presented through a forest plot, a prognostic risk model diagram, and a variable trajectory diagram. Multivariate Cox regression analysis incorporating the LASSO-derived risk score and clinical variables was performed to identify independent prognostic factors, with P<0.05 considered statistically significant.

The calculation for risk score was:

RiskScore=∑lCoefficient(genei)∗mRNAExpresion(genei)[1]

The R package ggplot2 was employed to illustrate the risk factors according to the LASSO risk score. Samples of HNSCC from the TCGA-HNSCC dataset were categorized into high- and low-risk groups according to the median LASSO risk score. The high-risk group consisted of samples with risk scores higher than the median, whereas the low-risk group consisted of samples with scores lower than the median. Kaplan-Meier (KM) curves (35) were generated to compare overall survival (OS) between the high- and low-risk groups. Cox proportional hazards regression analysis was performed to evaluate the prognostic significance of variables, and the results were presented as hazard ratios (HRs) with corresponding 95% confidence intervals (CIs). We estimated the 1-, 3-, and 5-year survival probability by constructing time-dependent receiver operating characteristic (ROC) curves and calculating the area under the curve (AUC) using the R package survivalROC (36) (Version 1.0.3.1).

Analysis of the prognostic risk model for HNSCC

The association among the risk score and clinical characteristics was demonstrated by a forest plot of the multivariate Cox regression findings. To illustrate the relationship among the risk score and clinical features, a nomogram (37) was created utilizing the R program rms (Version 6.7-1).

Making calibration curves that contrasted the expected and actual probabilities allowed us to assess the model’s prediction accuracy. Calibration analysis was conducted to assess the model’s reliability, whereas decision curve analysis (DCA) was executed using the R package ggDCA (38) (Version 1.1) to evaluate the clinical utility and discriminative capability of the LASSO risk score-based prognostic model for HNSCC.

Verification of the differential expression of key genes and ROC curve analysis

With the help of a group comparison plot, we were able to look at the TCGA-HNSCC and combined GEO datasets for genes that were differentially expressed across the control and HNSCC groups. The R package pROC (version 1.18.5) was utilized to construct the ROC curve and compute the AUC for evaluating the diagnostic significance of key gene expression. The AUC varies from 0.5 to 1, where higher values signify improved diagnostic performance (0.5–0.7: low, 0.7–0.9: moderate, >0.9: high).

Gene set variation analysis (GSVA)

Pathway enrichment variations among the samples were investigated using the R package GSVA (39) (Version 1.50.0). The analysis drew on the h.all.v7.4.symbols.gmt gene set derived from the MSigDB (32). Evaluated were functional enrichment variations between the high- and low-risk groups using P<0.05 as the significance level. Differential analysis of GSVA scores was performed using the limma package, and P values were adjusted for multiple testing using the BH method.

Calculation of the glycosylation and mitophagy score (G&Mscore)

The G&MScore for all samples was computed using the R package GSVA for ssGSEA (40), with the use of the expression matrices of important genes from the TCGA-HNSCC dataset. HNSCC samples were categorized into high- and low-score groups according to the median G&MScore. A comparative plot was created to examine expression differences between groups. Furthermore, KM curves (35) were plotted, and the R package survivalROC was utilized to produce time-dependent ROC curves and calculate the AUC for G&MScore and OS. The objective was to predict 1-, 3-, and 5-year survival probabilities for the HNSCC cohort. The AUC ranges from 0.5 to 1, with higher values indicating better diagnostic performance (0.5–0.7: low, 0.7–0.9: moderate, >0.9: high).

Analysis of immune-cell infiltration in the high- and low-risk groups

Activated CD8+ T cells, activated dendritic cells, gamma-delta T cells, natural killer cells, and regulatory T cells were among the immune cell types whose infiltration was evaluated using single-sample gene set enrichment analysis (ssGSEA) (40). Immune cell infiltration in the TCGA-HNSCC cohort was quantified using the ssGSEA algorithm based on the gene expression matrix. The relative abundance of 28 immune cell types was estimated for each sample. Differences in immune cell infiltration between the high- and low-risk groups were visualized using the ggplot2 package in R. For additional examination, immune cells exhibiting notable group differences were chosen. Spearman’s technique was used to evaluate the correlations between the important genes and immune cells, and the results were visualized using a correlation bubble plot.

Immunophenoscore (IPS) analysis

Immunogenicity denotes the capacity of an antigen or its epitopes to engage antigen-recognition receptors on T and B lymphocytes, thereby triggering humoral and/or cell-mediated immune responses. It is regulated by multiple genetic factors, including effector cell-related genes, major histocompatibility complex (MHC) molecules, immunomodulatory mediators, and immunosuppressive cell populations. Machine learning approaches can be applied to assess and quantify immunogenicity. The TCIA database (41) provides valuable resources for predicting CTLA-4 and PD-1 responsiveness. Accordingly, IPS data for HNSCC cases were retrieved from the TCGA-HNSCC dataset via the TCIA platform. Group-wise comparisons of IPS values between high- and low-risk cohorts were visualized using the R package ggplot2.

Statistical analysis

The R program, specifically version 4.2.2, was used for processing and analyzing the data. Continuous variables were assessed for normality prior to statistical analysis. If the continuous variable’s data followed a normal distribution, the Student’s t-test was utilized; otherwise, non-normally distributed data was evaluated using the Mann-Whitney U test. A correlation study called Spearman’s was run on the molecules to find their correlation coefficients. We fixed the significance level at P<0.05 and used two-sided P values.

Results

Technology roadmap

The overall analytical workflow of this study is shown in Figure 1. Briefly, transcriptomic data of HNSCC were obtained from the TCGA and GEO databases, and G&MRDEGs were identified. Functional enrichment analysis and GSEA were subsequently performed to explore the potential biological functions of these genes. Based on the selected candidate genes, a prognostic risk model was constructed and validated using survival analysis, ROC curves, and nomogram analysis. Finally, immune infiltration analysis and IPS evaluation were conducted to assess the potential clinical relevance of the prognostic model.

Figure 1 Flowchart for the comprehensive analysis. CNV, copy number variations; DCA, decision curve analysis; DEGs, differentially expressed genes; GO, Gene Ontology; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; HNSCC, head and neck squamous cell carcinomas; IPS, Immunophenoscore; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; G&MRDEGs, glycosylation- and mitophagy-related differentially expressed genes; G&MRGs, glycosylation-and mitophagy-related genes; G&MScore, Glycosylation and Mitophagy Score; ROC, receiver operating characteristic; SM, somatic mutation; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

Merging the HNSCC datasets

Batch effects between the GEO datasets (GSE107591 and GSE6631) were effectively corrected prior to data integration. As shown in Figure 2A,2B, the distribution of gene expression values across samples became more consistent after batch correction. In addition, PCA (23) demonstrated that samples from different datasets were more uniformly distributed following correction compared with the pre-correction state (Figure 2C,2D). These results indicate that batch-effect removal improved the comparability of the datasets and enabled the construction of a reliable integrated GEO dataset for subsequent analyses.

Figure 2 Removal of batch effects in the GSE107591 and GSE6631 datasets. (A,B) Boxplots showing the distribution of gene expression in the combined GEO datasets before (A) and after (B) correction of batch effects. (C,D) PCA plots illustrating sample clustering before (C) and after correction (D). GEO, Gene Expression Omnibus; PCA, principal component analysis.

G&MRDEGs within HNSCC

In TCGA-HNSCC dataset, 1,516 DEGs were found (|logFC| >2, adj. P<0.05), with 590 genes showing upregulation and 926 genes showing downregulation (Figure 3A). Furthermore, within the integrated GEO dataset, 2,089 DEGs were discovered, with 1,048 genes showing upregulation and 1,041 genes showing downregulation (|logFC| >0.25, P<0.05).

Figure 3 Differential gene expression analysis. (A) Volcano plot of the analysis of DEGs in HNSCC and control samples from the TCGA-HNSCC dataset. (B) DEGs in the TCGA-HNSCC dataset and combined GEO datasets, and the Venn diagram of genes related to G&MRGs. (C) Heat map of G&MRDEGs in the TCGA-HNSCC dataset. (D) Chromosomal localization map of G&MRDEGs. DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; G&MRDEGs, glycosylation- and mitophagy-related differentially expressed genes; G&MRGs, glycosylation- and mitophagy-related genes; HNSCC, head and neck squamous cell carcinoma; TCGA, The Cancer Genome Atlas.

To identify G&MRDEGs among the HNSCC samples, DEGs were first identified in the TCGA-HNSCC dataset (|logFC| >2, adj. P<0.05) and the combined GEO datasets (|logFC| >0.25, P<0.05). The intersection with G&MRGs yielded 23 G&MRDEGs (Table S1), which were visualized using a Venn diagram (Figure 3B). Using a heatmap, we compared the expression of these 23 G&MRDEGs in the TCGA-HNSCC dataset and showed the results in Figure 3C. Chromosomal locations were determined using RCircos. The evaluation revealed that many G&MRDEGs mapped to chromosomes 11 and 17, including MMPPARG3, MMP3, and SERPINH1 (Chr11) and BIRC5, COL1A1, and MAPT (Chr17) (Figure 3D).

CNVs and SMs of G&MRGs

To examine SMs of 866 G&MRGs in HNSCC samples from the TCGA-HNSCC dataset, conclusions of gene mutation investigation were aggregated and presented with the R package maftools (Figure 4A). The investigation identified missense mutations were the extremely prevalent of nine main types of SMs. Furthermore, single nucleotide polymorphism was the predominant mutation type; particularly prevalent single nucleotide variant in the HNSCC samples was C to T alterations.

Figure 4 CNV and somatic mutation analysis. (A) Presentation of SM of G&MRGs in HNSCC samples from the TCGA-HNSCC dataset. (B) Presentation of SM of G&MRDEGs in HNSCC samples from the TCGA-HNSCC dataset. (C,D) The CNV of G&MRDEGs are shown in the HNSCC samples from the TCGA-HNSCC dataset. CNV, copy number variation; G&MRDEGs, glycosylation- and mitophagy-related differentially expressed genes; G&MRGs, glycosylation- and mitophagy-related genes; HNSCC, head and neck squamous cell carcinoma; SM, somatic mutation; TCGA, The Cancer Genome Atlas.

The SM status of the 23 previously identified G&MRDEGs was analyzed in the HNSCC cohort, and the genes were ranked according to their mutation frequencies. As shown in Figure 4B, titin (TTN) exhibited the highest mutation frequency (36%).

To examine the CNVs of G&MRDEGs in HNSCC specimens derived from the TCGA-HNSCC dataset, CNV data were acquired, consolidated, and analyzed with GISTIC2.0 (26). A total of 23 G&MRDEGs with CNVs were identified in the HNSCC samples, and the detailed CNV status of these genes is provided in Table S2, while the corresponding mutation statuses were visualized (Figure 4C,4D).

GO and KEGG pathway enrichment in HNSCC

GO (27) enrichment analysis revealed that the G&MRDEGs were mainly involved in biological processes such as response to peptide, wound healing, and regulation of peptidase activity (Table 3), indicating potential roles in inflammatory responses and tissue remodeling (Figure 5A,5B). For cellular components, these genes were significantly enriched in the collagen-containing extracellular matrix and membrane raft, suggesting their association with tumor microenvironment-related structures (Figure 5C). For molecular function, the G&MRDEGs were significantly enriched in collagen binding, peptidase regulator activity, chaperone binding, protein C-terminus binding, and protease binding (Figure 5D), suggesting their involvement in extracellular matrix regulation, protease activity modulation, and protein interaction processes. KEGG pathway (28) analysis further demonstrated significant enrichment in the ECM-receptor interaction pathway (Figure 5E), implying that these genes may participate in HNSCC progression through extracellular matrix-mediated signaling pathways.

Figure 5 GO and KEGG enrichment analyses of G&MRDEGs. (A) Bubble plot of the results of GO and KEGG pathway enrichment analyses of G&MRDEGs, including BP, CC, MF, and biological pathways (KEGG). The abscissa shows the GO and KEGG terms. (B-D) Network diagrams of the results of GO enrichment analysis of G&MRDEGs. BP, biological process; CC, cellular component; G&MRDEGs, glycosylation- and mitophagy-related differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Gene set enrichment in HNSCC

We used GSEA to look at the TCGA-HNSCC dataset for expression profiles together with relevant biological processes, cellular components, and molecular functions to determine how gene expression levels affect HNSCC (Figure 6A).The research results (Table 4) demonstrated that the TCGA-HNSCC dataset has a substantial increase in processes and pathways including the assembly of collagen fibrils and other multimeric structures (Figure 6B), DNA damage, telomere stress-induced senescence (Figure 6C), HDAC deacetylation of histones (Figure 6D), and MET promotion of cell motility (Figure 6E).

Figure 6 GSEA for TCGA-HNSCC. (A) Mountain plot displaying four biological functions from the GSEA of the TCGA-HNSCC dataset. (B-E) GSEA results showing significant enrichment in the assembly of collagen fibrils and other multimeric structures (B), DNA damage, telomere stress-induced senescence (C), HDAC deacetylation of histones (D), and MET promotion of cell motility (E). GSEA, gene set enrichment analysis; HNSCC, head and neck squamous cell carcinoma; NES, normalized enrichment score; TCGA, The Cancer Genome Atlas.

The HNSCC prognostic risk model’s results

Clinical data and G&MRDEGs from the TCGA-HNSCC database were used in the development of the HNSCC predictive risk model by univariate Cox regression analysis. After using LASSO regression analysis with variables with P values <0.10, a six-gene model (PPARG, PIP, BIRC5, AR, ITGA5, ENO2) was identified (Figure 7A-7C). The following equation was used to obtain the risk score:

RiskScore = PPARG×0.2804−PIP×0.0882 + BIRC5×0.1435−AR×0.2824+ ITGA5×0.2063+ENO2×0.0249[2]

Figure 7 LASSO and Cox regression analyses. (A) Forest plot of G&MRDEGs based on univariate Cox regression analysis of prognosis. (B,C) Prognostic risk model plot (B) and variable trajectory plot (C) from the LASSO regression model. (D) Risk factor map for the prognostic LASSO model of key genes. (E) Prognostic KM curve comparing the high- and low-risk score groups from the LASSO model and OS in the HNSCC group. HR with 95% CI is shown. (F) Time-dependent ROC curve for the HNSCC group in the TCGA-HNSCC dataset. CI, confidence interval; HNSCC, head and neck squamous cell carcinoma; HR, hazard ratio; KM, Kaplan-Meier; LASSO, least absolute shrinkage and selection operator; MRDEGs, glycosylation- and mitophagy-related differentially expressed genes; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

The LASSO risk score was used to generate risk factor plots, which were generated using the ggplot2 package (Figure 7D). The findings demonstrated an increased mortality rate in the high-risk group, characterized by elevated expression levels of key genes in the prognostic model. A KM survival analysis, stratified by median OS, was performed utilizing LASSO risk scores and OS information obtained from the TCGA-HNSCC dataset (Figure 7E). The analysis indicated an important distinction in OS among high- and low-risk groups (P<0.001). Additionally, the time-dependent ROC curves (Figure 7F) exhibited moderate predictive accuracy (0.5< AUC <0.7) at 1-, 3-, and 5-year intervals.

After LASSO regression analysis, a multivariate Cox regression research was carried out to investigate how LASSO risk scores correlated with possible clinical outcomes for prognosis. An illustration of the findings was provided in the form of a forest plot (Figure 8A). To prove that the HNSCC risk model is useful for prognosis, we built a nomogram that shows the connections between genes and the disease by combining the LASSO risk score with clinical variables (Figure 8B). Among the factors tested, the LASSO risk score demonstrated the highest predictive value, whereas age was shown to be the least useful.

Figure 8 Prognostic analysis. (A) Forest plot of the univariate Cox regression analysis. (B) Nomogram constructed based on the multivariate Cox regression model. (C-E) Calibration curves of the prognostic risk model for HNSCC at 1 year (C), 3 years (D), and 5 years (E). (F-H) DCA plots of the prognostic risk model for HNSCC at 1 year (F), 3 years (G), and 5 years (H). DCA, decision curve analysis; HNSCC, head and neck squamous cell carcinoma.

The prognostic model underwent calibration at 1 year (Figure 8C), 3 years (Figure 8D), and 5 years (Figure 8E), with calibration curves illustrated (Figure 8C-8E). The findings indicated that the model exhibited optimal predictive performance at the 5-year mark.

We utilized DCA to assess the prognostic model’s treatment utility for HNSCC at 1, 3, and 5 years (Figure 8F-8H, respectively). Greater overall benefit is likely indicated by the model curve’s stability above the reference lines of “All positive” and “All negative” within a certain range. The prognostic model emerged victorious with the best clinical predictive value at 5 years, surpassing values at 3 and 1 years.

Key gene expression variations in HNSCC

We utilized the TCGA-HNSCC dataset to look at how key genes are expressed differently in various groups. The expression levels of six important genes (PPARG, PIP, BIRC5, AR, ITGA5, ENO2) differed significantly between the HNSCC and control samples (P<0.001), as shown in the results (Figure 9A). To create ROC curves using the TCGA-HNSCC dataset’s gene expression levels, the pROC package (36) in R was utilized (Figure 9B,9C). The results showed that when it came to differentiating HNSCC from control samples, BIRC5, ITGA5, and ENO2 had excellent diagnostic accuracy (AUC >0.9), while PPARG, PIP, and AR had moderate accuracy (0.7< AUC <0.9).

Figure 9 Validation of differential expression and ROC curve analysis of key genes. (A) Group comparison diagram of key genes in the HNSCC and control samples in the TCGA-HNSCC dataset. (B,C) ROC curves for the key genes PPARG, PIP, and BIRC5 (B) and AR, ITGA5, and ENO2 (C) in the TCGA-HNSCC dataset. (D) Group comparison plot of key genes in HNSCC and control samples in the combined GEO dataset. (E,F) ROC curves for the key genes PPARG, PIP, and BIRC5 (E) and AR, ITGA5, and ENO2 (F) in the combined GEO Datasets. Group comparisons were performed using an independent Student’s t-test. Statistical significance was indicated as ***, P<0.001. GEO, Gene Expression Omnibus; HNSCC, head and neck squamous cell carcinoma; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

To assess the variations in gene expression levels across the GEO datasets, we used a differential expression analysis. The findings (Figure 9D) indicated that five pivotal genes (PPARG, BIRC5, AR, ITGA5, and ENO2) were markedly elevated in HNSCC samples relative to control samples (P<0.001). ROC curves were produced utilizing the pROC package in R to assess the diagnostic capability of these genes (Figure 9E,9F). The findings demonstrated that PPARG, BIRC5, AR, ITGA5, and ENO2 displayed intermediate diagnostic accuracy (0.7< AUC <0.9), while PIP exhibited lesser accuracy (0.5< AUC <0.7).

Gene set variation in HNSCC

The h.all.v7.4.symbols.gmt gene set was analyzed using GSVA (39) to compare the TCGA-HNSCC dataset’s high-risk and low-risk groups (Table 5). The ten pathways exhibiting the greatest logFC values with positive enrichment (P<0.05) and the ten pathways demonstrating negative enrichment were found. A heatmap (Figure 10A) was used to compare the expression variations of these 20 pathways in both groups at risk.

Figure 10 GSVA analysis. Heat map (A) and group comparison map (B) of GSVA results in the high-risk and low-risk groups in the TCGA-HNSCC dataset. Group comparisons were performed using an independent Student’s t-test. Statistical significance was indicated as *, P<0.05; **, P<0.01; ***, P<0.001. GSVA, gene set variation analysis; HNSCC, head and neck squamous cell carcinoma; TCGA, The Cancer Genome Atlas.

We utilized the Mann-Whitney U test to look more closely at the differences that were found to see if they were statistically significant. This results were then displayed in a group comparison diagram (Figure 10B). GSVA revealed statistically significant differences (P<0.05) in multiple pathways, including epithelial-mesenchymal transition, angiogenesis, glycolysis, hypoxia, mechanistic target of rapamycin complex 1 (mTORC1) signaling, unfolded protein response, MYC targets v2, G2/M checkpoint, E2F targets, MYC targets v1, pancreatic beta-cell function, spermatogenesis, KRAS signaling upregulation, interleukin-2/signal transducer and activator of transcription 5 (IL2-STAT5) signaling, apoptosis, peroxisome activity, allograft rejection, fatty acid metabolism, bile acid metabolism, and KRAS signaling downregulation.

Construction of the G&MScore

The G&MScore for all samples, derived from the expression of six pivotal genes in the TCGA-HNSCC database, was computed utilizing the ssGSEA (40) algorithm. The disparity in G&MScore among HNSCC and control samples was illustrated by a group comparison diagram (Figure 11A), revealing a statistically significant difference (P<0.05).

Figure 11 ssGSEA score analysis. (A) Group comparison plot of G&MScores in the HNSCC and control groups in the TCGA-HNSCC dataset. (B) Group comparison diagram of key genes in the high- and low-score groups of the HNSCC samples. (C) The prognostic KM curve illustrates the differences in OS between the high and low G&MScore groups within the HNSCC cohort. (D) Time-dependent ROC curve analysis of the HNSCC group using data from the TCGA-HNSCC dataset. Group comparisons were performed using an independent Student’s t-test. Statistical significance was indicated as *, P<0.05; **, P<0.01; ***, P<0.001. Hazard ratios (HR) with corresponding 95% CI are shown. CI, confidence interval; HNSCC, head and neck squamous cell carcinoma; KM, Kaplan-Meier; G&MScore, glycosylation and mitophagy score; OS, overall survival; ROC, receiver operating characteristic; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

The HNSCC samples were split into two groups based on the median G&MScore: high-score and low-score. These groups had distinct expression levels for six key genes (Figure 11B). We found that the expressions of all six genes differed significantly (P<0.05). ENO2, PPARG, PIP, ITGA5, and AR exhibited extremely significant variations among the groups (P<0.001), but BIRC5 was also significant (P<0.01).

The predictive value of the G&MScore for overall survival in the HNSCC group was evaluated using a KM survival analysis (Figure 11C). The results indicated a significant difference in OS rates between the groups with high and low scores (P<0.05). Furthermore, the ROC curves produced for 1-, 3-, and 5-year survival projections (Figure 11D) demonstrated that the HNSCC phenotypic scoring model displayed low predictive accuracy (0.5< AUC <0.7).

Infiltration of immune cells in high- and low-risk categories

A group comparison plot was utilized to illustrate the variation in the expression of invading immune cells across distinct groups. It was found that 14 types of immune cells (P<0.05) were statistically different between the groups (Figure 12A). These cells included activated B cells, activated CD8⁺ T cells, effector memory CD8⁺ T cells, immature B cells, memory B cells, type 1 T helper cells, type 17 T helper cells, CD56 bright natural killer cells, eosinophils, mast cells, myeloid-derived suppressor cells, monocytes, and neutrophils. Correlation bubble plots showed the relationships between important genes and immune-cell infiltration (Figure 12B,12C). The correlation bubble plot indicated that the principal genes displayed a robust positive connection with immune cells in both the low- and high-risk HNSCC cohorts. ITGA5 demonstrated the most robust positive connection with central memory CD8⁺ T cells in the low-risk cohort (r=0.610, P<0.05), while AR displayed the strongest positive correlation with central memory CD8⁺ T cells in the high-risk cohort (r=0.468, P<0.05).

Figure 12 Analysis of immune-cell infiltration in the risk groups using the ssGSEA algorithm. (A) Comparison of immune cells between the low-risk and high-risk HNSCC groups from the TCGA-HNSCC dataset. The bubble plot shows the abundance of immune-cell infiltration and the correlations between key genes in the low- (B) and high-risk (C) HNSCC groups. ns, no significance; *, P<0.05; **, P<0.01; ***, P<0.001. HNSCC, head and neck squamous cell carcinoma; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

Association between the risk score and IPS

We acquired IPS data pertaining to HNSCC samples from the TCIA database to evaluate the possible influence of immunotherapy on risk groups. We utilized the ggplot2 tool in R to illustrate group comparisons of the IPS among the two groups of LASSO risk score cohorts in the HNSCC data (Figure 13A-13D). The results indicated an important distinction in both positive and negative IPS for PD1 and CTLA4 among the two risk categories (P<0.05).

Figure 13 IPS analysis. (A-D) Group comparison of the IPS between the high-risk and low-risk groups in the HNSCC samples from the TCGA-HNSCC dataset. *, P<0.05; **, P<0.01; ***, P<0.001. HNSCC, head and neck squamous cell carcinoma; IPS, Immunophenoscore; TCGA, The Cancer Genome Atlas.

Discussion

HNSCC is a very aggressive malignant neoplasm that significantly impairs quality of life and is linked to unfavorable survival outcomes. Notwithstanding comprehensive research and advancements in treatment modalities during recent decades, the prognosis for HNSCC remains inadequate, and its incidence continues to rise (1,42). Using bioinformatics techniques, this work dug deep into the expression patterns and therapeutic relevance of G&MRGs in HNSCC. The HNSCC and control samples showed a notable variation in the expression of six genes (PPARG, PIP, BIRC5, AR, ITGA5, and ENO2) with a P value significantly lower than 0.001. The roles of these genes in cancer were further validated, revealing their potential roles in prognosis prediction and personalized treatment.

The findings of this investigation demonstrated that the transcript levels of BIRC5, ITGA5, and ENO2 were particularly effective in differentiating HNSCC samples from control samples. Notably, these genes play important roles in tumor development. BIRC5 has been recognized as an anti-apoptotic protein in multiple malignancies, with elevated expression correlating with tumor cell viability and resistance to treatment (43,44). ITGA5 plays a crucial role in ECM-receptor interactions, and its high expression may promote tumor cell invasion and metastasis (45). ENO2, a glycolytic enzyme, is related to metabolic reprogramming in tumor cells (46). Furthermore, PPARG, PIP, and AR demonstrated moderate classification accuracy in the present study (0.7< AUC <0.9). These results demonstrate the importance of these genes in the biological processes that cause HNSCC and suggest a possible relationship between their expression levels and the development of the disease.

The predictive risk model, based on six key genes, demonstrated patients classified as low-risk had significantly improved survival rates over the ones categorized as the high-risk category, a finding supported by validation results from the GEO cohort. Thus, these gene expression levels may be useful indicators for HNSCC prognosis, allowing for more accurate survival prediction and the creation of tailored therapy regimens.

The enrichment analysis performed in the present investigation indicated that the ECM–receptor interaction pathway was considerably enriched in HNSCC. This finding aligns with the invasive and metastatic characteristics of HNSCC. ECM–receptor interactions not only play key roles in cell adhesion and migration but also promote tumor cell growth and survival by activating signaling pathways (47). This finding offers significant insights into the cellular processes of HNSCC and establishes a theoretical foundation for the development of treatments aimed at ECM–receptor interactions. In the group with low risk, ITGA5 was the most significantly connected variable with CD8⁺ T cells in central memory, according to immune-cell infiltration analysis. The highest level of correlation was seen between AR and the high-risk category. ITGA5 is involved in cell adhesion and signaling processes and can promote the recruitment and retention of central memory CD8⁺ T cells within the tumor microenvironment, thereby enhancing their antitumor activity. This increased infiltration can lead to improved immune surveillance and destruction of tumor cells, thereby improving prognosis (48). Moreover, ITGA5 interacts with fibronectin in the ECM, influencing the structure and composition of the ECM. These changes can create a more favorable environment for immune-cell infiltration and function, further boosting anti-tumor responses (49). Directly influencing the expression of genes that regulate the immune system, AR can impact immune cells through both DNA binding-dependent and -non-dependent mechanisms (50). Thus, it appears that the immunological microenvironment can be modulated by the activity of key genes, which in turn affects the outcome of HNSCC. For example, ITGA5 may enhance antitumor immune responses by increasing the infiltration of central memory CD8⁺ T cells, thereby improving prognosis. In contrast, AR may affect immune cell function and distribution through different mechanisms that influence tumor progression.

Beyond the analyses of gene expression and functional enrichment, we further characterized the SM and CNV profiles of G&MRGs. Missense mutations were identified as the predominant mutation type, while single-nucleotide variations were mainly characterized by C>T substitutions, suggesting a potentially distinct mutational pattern in HNSCC, consistent with previous genomic studies (51). In addition, several candidate genes exhibited varying degrees of copy number amplification or deletion, reflecting a certain degree of genomic instability in HNSCC (52). Interestingly, a subset of G&MRDEGs appeared to cluster on chromosomes 11 and 17, regions that have been reported to undergo recurrent structural alterations in multiple malignancies (53). Collectively, these observations provide complementary genomic evidence supporting the potential involvement of G&MRGs in HNSCC. Nevertheless, the underlying molecular mechanisms remain to be elucidated and require further experimental validation.

The systematic bioinformatics analysis conducted in this study revealed the clinical significance of G&MRGs in HNSCC. To our knowledge, this study is the first in which conclusions were derived from analyzing the molecular mechanisms and clinical relevance of HNSCC from multiple perspectives and levels, providing a new theoretical basis and identifying key genes for the diagnosis, treatment, and prognostic evaluation of this disease. Nonetheless, this study possesses certain shortcomings that require attention. First, the results were primarily obtained using public database data, which were not subjected to further experimental validation. Second, the sample size was relatively small. Therefore, future large-scale studies are required to verify these findings. Moreover, forthcoming research should concentrate on integrating experimental biology techniques to investigate the precise processes by which these key genes function in HNSCC and assess their viability as potential treatments.

In conclusion, we meticulously examined the critical functions of glycosylation and mitophagy in HNSCC, offering novel insights and a theoretical framework for comprehending the molecular causes of HNSCC and enhancing the formulation of personalized treatment options for HNSCC. This study’s results provide a foundation for future experimental investigations to develop tumor prognosis biomarkers that could help doctors diagnose, predict, and treat HNSCC patients more effectively.

Conclusions

We identified and validated a prognostic gene signature associated with HNSCC and developed a predictive model with potential clinical utility. These findings may provide insights into prognostic stratification and individualized management of HNSCC. Further external and experimental validation is warranted.

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-1-0011/rc

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

Funding: This study was supported by the Zhejiang Provincial Medical and Health Science and Technology Plan Project (No. 2025KY1948).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0011/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. Chow LQM. Head and Neck Cancer. N Engl J Med 2020;382:60-72. [Crossref] [PubMed]
2. Goel B, Tiwari AK, Pandey RK, et al. Therapeutic approaches for the treatment of head and neck squamous cell carcinoma-An update on clinical trials. Transl Oncol 2022;21:101426. [Crossref] [PubMed]
3. Chen H, Zhang Y, Shen Y, et al. Deficiency of N-linked glycosylation impairs immune function of B7-H6. Front Immunol 2023;14:1255667. [Crossref] [PubMed]
4. Sharma P, Zhang X, Ly K, et al. Hyperglycosylation of prosaposin in tumor dendritic cells drives immune escape. Science 2024;383:190-200. [Crossref] [PubMed]
5. Chandler KB, Costello CE, Rahimi N. Glycosylation in the Tumor Microenvironment: Implications for Tumor Angiogenesis and Metastasis. Cells 2019;8:544. [Crossref] [PubMed]
6. Shin YY, Seo Y, Oh SJ, et al. Melatonin and verteporfin synergistically suppress the growth and stemness of head and neck squamous cell carcinoma through the regulation of mitochondrial dynamics. J Pineal Res 2022;72:e12779. [Crossref] [PubMed]
7. Bertola N, Degan P, Cappelli E, et al. Mutated FANCA Gene Role in the Modulation of Energy Metabolism and Mitochondrial Dynamics in Head and Neck Squamous Cell Carcinoma. Cells 2022;11:2353. [Crossref] [PubMed]
8. Liao C, An J, Tan Z, et al. Changes in Protein Glycosylation in Head and Neck Squamous Cell Carcinoma. J Cancer 2021;12:1455-66. [Crossref] [PubMed]
9. Guo LL, Yang XH, Jin HZ, et al. TFAP2A induces cisplatin resistance via BNIP3-mediated mitophagy in non-small cell lung cancer. Transl Cancer Res 2025;14:5353-69. [Crossref] [PubMed]
10. Yanan L, Hui L, Zhuo C, et al. Comprehensive analysis of mitophagy in HPV-related head and neck squamous cell carcinoma. Sci Rep 2023;13:7480. [Crossref] [PubMed]
11. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
12. Goldman MJ, Craft B, Hastie M, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol 2020;38:675-8. [Crossref] [PubMed]
13. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 2007;23:1846-7. [Crossref] [PubMed]
14. Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res 2013;41:D991-5. [Crossref] [PubMed]
15. Verduci L, Ferraiuolo M, Sacconi A, et al. The oncogenic role of circPVT1 in head and neck squamous cell carcinoma is mediated through the mutant p53/YAP/TEAD transcription-competent complex. Genome Biol 2017;18:237. [Crossref] [PubMed]
16. Sacconi A, Donzelli S, Pulito C, et al. TMPRSS2, a SARS-CoV-2 internalization protease is downregulated in head and neck cancer patients. J Exp Clin Cancer Res 2020;39:200. [Crossref] [PubMed]
17. Kuriakose MA, Chen WT, He ZM, et al. Selection and validation of differentially expressed genes in head and neck cancer. Cell Mol Life Sci 2004;61:1372-83. [Crossref] [PubMed]
18. Leek JT, Johnson WE, Parker HS, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012;28:882-3. [Crossref] [PubMed]
19. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
20. Stelzer G, Rosen N, Plaschkes I, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics 2016;54:1.30.1-1.30.33.
21. Qiu P, Guo Q, Yao Q, et al. Characterization of Exosome-Related Gene Risk Model to Evaluate the Tumor Immune Microenvironment and Predict Prognosis in Triple-Negative Breast Cancer. Front Immunol 2021;12:736030. [Crossref] [PubMed]
22. Bi J, Bi F, Pan X, et al. Establishment of a novel glycolysis-related prognostic gene signature for ovarian cancer and its relationships with immune infiltration of the tumor microenvironment. J Transl Med 2021;19:382. [Crossref] [PubMed]
23. Ben Salem K, Ben Abdelaziz A. Principal Component Analysis (PCA). Tunis Med 2021;99:383-9.
24. Zhang H, Meltzer P, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinformatics 2013;14:244. [Crossref] [PubMed]
25. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
26. Li J, Miao B, Wang S, et al. Hiplot: a comprehensive and easy-to-use web service for boosting publication-ready biomedical data visualization. Brief Bioinform 2022;23:bbac261. [Crossref] [PubMed]
27. Mi H, Muruganujan A, Ebert D, et al. PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res 2019;47:D419-26. [Crossref] [PubMed]
28. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28:27-30. [Crossref] [PubMed]
29. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
30. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
31. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics 2013;29:1830-1. [Crossref] [PubMed]
32. Liberzon A, Subramanian A, Pinchback R, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 2011;27:1739-40. [Crossref] [PubMed]
33. T, T. A Package for Survival Analysis in R. R package version 3.4-0 (2022). Available online: https://cran.r-project.org/package=survival
34. Engebretsen S, Bohlin J. Statistical predictions with glmnet. Clin Epigenetics 2019;11:123. [Crossref] [PubMed]
35. Rich JT, Neely JG, Paniello RC, et al. A practical guide to understanding Kaplan-Meier curves. Otolaryngol Head Neck Surg 2010;143:331-6. [Crossref] [PubMed]
36. Park SH, Goo JM, Jo CH. Receiver operating characteristic (ROC) curve: practical review for radiologists. Korean J Radiol 2004;5:11-8. [Crossref] [PubMed]
37. Van Calster B, Wynants L, Verbeek JFM, et al. Reporting and Interpreting Decision Curve Analysis: A Guide for Investigators. Eur Urol 2018;74:796-804. [Crossref] [PubMed]
38. Wu J, Zhang H, Li L, et al. A nomogram for predicting overall survival in patients with low-grade endometrial stromal sarcoma: A population-based analysis. Cancer Commun (Lond) 2020;40:301-12. [Crossref] [PubMed]
39. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
40. Xiao B, Liu L, Li A, et al. Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma. Front Oncol 2020;10:607622. [Crossref] [PubMed]
41. Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
42. Mereiter S, Balmaña M, Campos D, et al. Glycosylation in the Era of Cancer-Targeted Therapy: Where Are We Heading? Cancer Cell 2019;36:6-16. [Crossref] [PubMed]
43. Li F, Aljahdali I, Ling X. Cancer therapeutics using survivin BIRC5 as a target: what can we do after over two decades of study? J Exp Clin Cancer Res 2019;38:368. [Crossref] [PubMed]
44. Liu SH, Hong Y, Markowiak S, et al. BIRC5 is a target for molecular imaging and detection of human pancreatic cancer. Cancer Lett 2019;457:10-9. [Crossref] [PubMed]
45. Kuninty PR, Bansal R, De Geus SWL, et al. ITGA5 inhibition in pancreatic stellate cells attenuates desmoplasia and potentiates efficacy of chemotherapy in pancreatic cancer. Sci Adv 2019;5:eaax2770. [Crossref] [PubMed]
46. Gao L, Yang F, Tang D, et al. Mediation of PKM2-dependent glycolytic and non-glycolytic pathways by ENO2 in head and neck cancer development. J Exp Clin Cancer Res 2023;42:1. [Crossref] [PubMed]
47. Eble JA, Niland S. The extracellular matrix in tumor progression and metastasis. Clin Exp Metastasis 2019;36:171-98. [Crossref] [PubMed]
48. Loo CS, Gatchalian J, Liang Y, et al. A Genome-wide CRISPR Screen Reveals a Role for the Non-canonical Nucleosome-Remodeling BAF Complex in Foxp3 Expression and Regulatory T Cell Function. Immunity 2020;53:143-157.e8. [Crossref] [PubMed]
49. Lu L, Gao Y, Huang D, et al. Targeting integrin alpha5 in fibroblasts potentiates colorectal cancer response to PD-L1 blockade by affecting extracellular-matrix deposition. J Immunother Cancer 2023;11.
50. Ben-Batalla I, Vargas-Delgado ME, von Amsberg G, et al. Influence of Androgens on Immunity to Self and Foreign: Effects on Immunity and Cancer. Front Immunol 2020;11:1184. [Crossref] [PubMed]
51. Cheng H, Jiang M, Kim S, et al. Uncovering chromatin factor landscapes in head and neck squamous cell carcinoma. Oral Oncol 2025;170:107687. [Crossref] [PubMed]
52. Dong Y, Qing M, Zhang Y, et al. Multi-omics analysis reveals the association between intratumoral bacteria and somatic mutational signatures in oral squamous cell carcinoma. J Transl Med 2025;24:23. [Crossref] [PubMed]
53. Vitale I, Cereda M, Galluzzi L. Epigenetic drivers of chromosomal instability. Trends Cancer 2026;12:95-8. [Crossref] [PubMed]