E2F expression profiling-based subtypes in head and neck squamous cell carcinoma: clinical relevance, prognostic implications, and personalized therapy

Abstract

Background

Head and neck squamous cell carcinoma (HNSCC) is a heterogeneous malignancy with poor prognosis. Dysregulation of E2F transcription factors (E2Fs), which control cell proliferation and apoptosis, is implicated in HNSCC pathogenesis. This study explores HNSCC molecular heterogeneity via E2Fs expression, identifies distinct subtypes, and develops a prognostic model that integrates gene expression, immune infiltration, and drug sensitivity.

Methods

We analyzed the TCGA-HNSC dataset (n = 494) and classified samples based on the expression of eight E2Fs using ConsensusClusterPlus. The optimal number of clusters (k = 2) was determined with the getOptK() function, which assesses cluster stability via internal consistency metrics. Differentially expressed genes between subtypes were identified with limma, and functional annotation was performed using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses. A prognostic model was constructed using LASSO regression on genes significant in univariate Cox analysis and validated in an independent GSE41613 dataset (n = 97). Immune cell infiltration was estimated using CIBERSORT, and drug sensitivity predicted via pRRophetic. Confounding factors such as HPV and smoking status were not included due to incomplete data.

Results

Two distinct E2F-based subtypes emerged. Cluster 1, characterized by lower E2Fs expression, exhibited poorer overall survival (log-rank, p = 0.035) and was enriched in genes related to epidermal development, keratinocyte differentiation, and IL-17 signaling. In contrast, Cluster 2 showed higher E2Fs expression, better survival, and enrichment in genes associated with DNA replication and repair. Notably, high-risk patients demonstrated increased infiltration of M0 and M2 macrophages (p < 0.05), suggesting an immunosuppressive tumor microenvironment that adversely affects prognosis. Our seven-gene prognostic model (AREG, CXCL14, FAM83E, FDCSP, ARHGAP4, EPHX3, and SPINK6) exhibited robust performance with AUCs of 0.692, 0.673, and 0.679 for 1-, 3-, and 5-year survival, a C-index of 0.66, and good calibration. High-risk patients also showed greater sensitivity to targeted agents such as pazopanib and imatinib.

Conclusions

These findings reveal two distinct E2F-based molecular subtypes of HNSCC that differ in prognosis, functional pathways, immune infiltration, and drug sensitivity. The prognostic model offers valuable risk stratification and identifies potential biomarkers and therapeutic targets, warranting further experimental and clinical validation.

Introduction

Head and neck squamous cell carcinoma (HNSCC) is a prevalent and aggressive cancer affecting the oral cavity, pharynx, and larynx, accounting for approximately 4.41% of global cancer cases [1]. Despite advancements in surgical, radiation, and chemotherapy treatments, the prognosis of HNSCC remains poor, with a five-year survival rate of less than 50% and a high recurrence rate [2, 3]. The complexity of HNSCC is compounded by its biological heterogeneity, driven by various etiological factors such as tobacco use, alcohol consumption, and human papillomavirus (HPV) infection [4].

This biological variability presents a significant challenge to effective treatment [5]. Additionally, intrinsic and acquired resistance mechanisms, driven by pathways such as EGFR signaling and immune evasion, often limit the success of current therapies. Although emerging treatments such as targeted therapies and immunotherapies show promise, their efficacy varies widely among patients, underscoring the urgent need for personalized treatment strategies. A more precise understanding of the molecular mechanisms, especially those driving resistance, is crucial for improving outcomes.

E2F transcription factors (E2Fs) are crucial regulators of the cell cycle, overseeing processes such as proliferation, differentiation, and apoptosis [6]. In normal cells, their activity is tightly controlled by regulators like the retinoblastoma protein. However, in cancer, genetic mutations, epigenetic alterations, or oncogenic interactions disrupt this regulation, leading to aberrant E2F activation that drives uncontrolled growth and tumor progression [7]. In HNSCC, frequent abnormalities in the E2F pathway contribute significantly to tumorigenesis by promoting the expression of genes involved in DNA replication, cell cycle progression, and apoptosis inhibition [8]. Recent studies have further shown that dysregulated E2F family members not only underlie tumor heterogeneity in HNSCC and other cancers, but also correlate with clinical features such as tumor stage, grade, and prognosis [9, 10]. Furthermore, E2F expression profiles correlate with immune cell infiltration in the tumor microenvironment (TME), highlighting their role in shaping both molecular and immune landscapes [11]. Despite this, few studies have systematically integrated E2F expression patterns with molecular subtyping, prognosis, immune features, and drug response in HNSCC.

Given the variability of E2F expression and its critical role in tumor biology, we hypothesize that distinct E2F expression profiles may define HNSCC subtypes with differing clinical outcomes and therapeutic responses. By elucidating the role of E2Fs in tumor progression and immune modulation, this molecular subtyping strategy could not only serve as a biomarker for prognosis and drug sensitivity but also pave the way for personalized treatment strategies and the discovery of novel therapeutic targets.

Materials and methods

Collection of data and preprocessing

The TCGA-HNSC dataset from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga) was obtained using the “TCGAbiolinks” R package (version 2.30.0) [12]. The dataset comprised 494 tumor samples with complete gene expression profiles data and overall survival (OS) details. Baseline clinical characteristics of patients in the TCGA-HNSC are presented in Supplementary Table S1. Due to the high proportion of missing data for HPV and smoking status in the TCGA-HNSC dataset, they were not incorporated into the prognostic nomogram, but were considered in subgroup analyses for prognostic evaluation. The study focused on eight E2Fs (E2F1 to E2F8), which play critical roles in cell cycle regulation and tumorigenesis.

For validation, we utilized the GSE41613 dataset [13], which comprises 97 oral squamous cell carcinoma patient samples retrieved from the Gene Expression Omnibus (GEO) database on the GPL570 platform. Data processing involved removing control probes, excluding probes that mapped to multiple genes, and calculating the median expression values for genes represented by multiple probes.

HNSCC sample subtyping and prognostic analysis

To identify robust molecular subtypes of HNSCC, we performed consensus clustering using the ConsensusClusterPlus [14] R package. The optimal number of clusters was determined with the getOptK() function, based on the expression profiles of eight E2Fs. Clustering was performed using the k-means algorithm with Euclidean distance, a widely used method for high-dimensional gene expression data [15]. The analysis included 1000 iterations with 80% item resampling to ensure stable subtypes. Kaplan–Meier survival curves were plotted, and survival differences between subtypes were assessed using the log-rank test.

Analysis of differential gene expression and functional enrichment

Differentially expressed genes (DEGs) analysis between HNSCC subtypes was performed using the limma [16] R package (v3.58.1). DEGs were identified based on|log2(FC)| > 1 and a q-value < 0.01, as determined by the Benjamini–Hochberg correction, to ensure the inclusion of genes with biologically meaningful expression differences. The ClusterProfiler [17] R package was utilized to perform enrichment analysis of DEGs for Gene Ontology (GO) [18] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [19] pathways. An initial p-value cutoff of 0.05 was applied, followed by multiple testing correction using the Benjamini–Hochberg method, with a final q-value threshold of 0.05. Results of the enrichment analysis were visualized using the ggplot2 R package.

Construction and validation of the prognostic model

The coxph function from the R survival [20] package was used to perform univariate Cox regression analysis on all DEGs. Genes with a p-value < 0.05 were identified as significant prognostic indicators for HNSCC. Given that Cox regression is designed to assess survival associations rather than expression differences, a fold-change threshold was not employed at this stage. The identified prognostic factors were visualized using a forest plot, which was generated using the forestplot R package. Subsequently, a prognostic model was constructed using least absolute shrinkage and selection operator (LASSO) regression via the glmnet [21] R package, incorporating significant prognostic factors identified from the univariate Cox analysis. A risk score for each patient was determined by calculating the weighted sum of gene expression levels, where 𝛽𝑖 denotes the weight coefficient and 𝜒𝑖 denotes the gene expression level.

$$\:Risk\:Score = \sum\limits_{i = 0}^n {\beta i*\chi i} $$

Patients with high or low scores were categorized on the basis of their median risk scores. Survival analysis utilized the survminer [22] R package, while receiver operating characteristic (ROC) curve analysis employed the timeROC [23] R package. Cox regression analyses, both univariate and multivariate, were conducted to assess if the risk score independently predicts prognosis. The same formula was used in the validation cohort to validate the model and calculate risk scores. A nomogram was created using the R package rms to assess the predictive accuracy of the model. Calibration curves were generated to evaluate the performance of the nomogram.

Comparison with traditional TNM staging system

To assess the prognostic performance of the constructed risk model in comparison with the traditional TNM staging system, clinical staging information based on the AJCC 7th edition was extracted for patients in the TCGA-HNSC cohort. Patients were categorized into early-stage (stage I/II) and advanced-stage (stage III/IV) groups. Kaplan–Meier survival analysis and log-rank tests were conducted to evaluate differences in OS between these groups.

Examination of immune cell infiltration

CIBERSORT classifies 22 immune cell types and has been widely validated across various tissue types, particularly in TME studies, where it correlates well with clinical outcomes. To explore the immune landscape of the E2F-based HNSCC subtypes, we used the CIBERSORT algorithm [24] to estimate the relative abundance of immune cells in the TME, implemented via the IOBR [25] R package. Spearman rank correlation analysis was conducted to examine the relationships between the risk score, its constituent genes, and immune cell proportions.

Analysis of drug sensitivity

We used the pRRophetic [26] R package to predict drug sensitivity in HNSCC patients, utilizing data from the Genomics of Drug Sensitivity in Cancer database [27] (http://www.cancerrxgene.org). The package employs a ridge regression model trained on gene expression and drug sensitivity data from GDSC to estimate the half-maximal inhibitory concentration (IC50) values for each drug in TCGA-HNSC samples. To reduce the false positive rate, the Benjamini–Hochberg method was used for multiple testing correction throughout the analysis. Differential expression analysis and comparisons of log-transformed IC50 values between high- and low-risk groups were both subjected to this correction. An adjusted p-value < 0.05 was considered statistically significant.

Statistical analysis

Statistical analyses were conducted using R software (version 4.2.0; R Foundation for Statistical Computing, Vienna, Austria). All results were considered statistically significant when the p-value was < 0.05. The notations used are: **** for p < 0.0001, *** for p < 0.001, ** for p < 0.01, * for p < 0.05, and ns for non-significant results.

Results

Identification and characterization of E2F-based subtypes in TCGA-HNSC

Consensus clustering of TCGA-HNSC samples was performed based on the expression of eight E2Fs. Based on Euclidean distance and k-means clustering, the optimal number of clusters was determined to be two, supported by the consensus matrix heatmap, resulting in the classification of samples into two molecular subtypes (Fig. 1A). Figure 1B illustrates a statistically significant difference in prognosis among the subtypes. A significantly poorer prognosis was observed in Cluster 1 compared to Cluster 2 (HR = 1.40; 95% CI: 1.02–1.90; p = 0.035). The expression levels of the eight E2Fs in the two subtypes, along with the clinical features associated with TCGA-HNSC samples, are shown in Fig. 1C. Expression levels of E2Fs in Cluster 1 were significantly lower than those in Cluster 2.

Fig. 1

Identification and characterization of E2F-based subtypes in TCGA-HNSC. (A) Two subtypes were identified. (B) Overall survival differed significantly between the two subtypes (p < 0.05). (C) E2F expression was significantly lower in Cluster 1 compared to Cluster 2. TCGA, The Cancer Genome Atlas; HNSC, head and neck squamous cell carcinoma

Analysis of differential gene expression and functional enrichment

Analysis of DEGs between the two clustered subtypes revealed 211 genes with differential expression, including 102 upregulated in Cluster 1 and 109 downregulated genes in Cluster 2. DEGs and principal component analysis (PCA) results for the two E2F-based subtypes are shown in Supplementary Fig. S1. Functional enrichment analysis of DEGs revealed a significant correlation with the intermediate filament cytoskeleton. The GO enrichment analysis initially identified several redundant biological processes. To improve clarity, we employed a semantic similarity approach using Revigo (http://revigo.irb.hr/) to consolidate similar terms and reduce redundancies, thereby ensuring clearer biological interpretations. Additionally, the GO enrichment analysis of upregulated genes in Cluster 1 revealed significant involvement in epidermal development, as well as enrichment in cellular components such as intermediate filaments and keratin envelopes. In terms of molecular function, these genes were associated with skin epidermal structural components and serine-type endopeptidase activity. KEGG pathway analysis further indicated significant enrichment in the IL-17 signaling pathway (Fig. 2A). Conversely, GO and KEGG analyses of upregulated genes in Cluster 2 showed that these genes are primarily involved in DNA replication, alterations in DNA topology, and double-strand break repair. These genes were enriched in cellular components such as nuclear chromosomal regions, while their molecular functions were primarily associated with single-stranded DNA binding, ATP-dependent activity on DNA, and helicase activity. (Fig. 2B). Additionally, KEGG pathway analysis revealed significant enrichment in the cell cycle and DNA replication pathways.

Fig. 2

Functional enrichment of genes differentially expressed between HNSCC subtypes. (A) Dot plot of GO and KEGG enrichment for genes upregulated in Cluster (1) (B) Dot plot of GO and KEGG enrichment for genes upregulated in Cluster (2) Dot size: gene count; HNSCC, Head and neck squamous cell carcinoma; GO, Gene ontology; BP, Biological process; CC, Cellular component; MF, Molecular function; KEGG, Kyoto encyclopedia of genes and genomes

Prognostic factor identification and prognostic model construction

Univariate Cox analysis identified 14 genes as significant prognostic factors (Fig. 3A). The optimal λ value resulted in the identification of seven genes (Fig. 3B-E) with corresponding coefficients. The risk score was determined for each patient by using these coefficients. Patients were divided into high- and low-score groups based on the median score. The survival probabilities differed significantly between high- and low-score groups (log-rank test, p < 0.001). Time-dependent ROC analysis was used to evaluate the predictive accuracy of the model for 1-, 3-, and 5-year survival rates. The areas under the curve (AUCs) at 1, 3, and 5 years were 0.692, 0.673, and 0.679, respectively (Fig. 4A-E), with a C-index of 0.66. The calibration curve (Fig. 4F) further demonstrated good agreement between the predicted and observed survival probabilities, indicating that the model predictions were well-calibrated in the training cohort.

Fig. 3

Prognostic factor identification and risk score model construction. (A) Univariate Cox regression analysis identifies 14 genes as significant prognostic factors. (B)-(E) Risk score formula constructed using LASSO regression. LASSO, least absolute shrinkage and selection operator

Fig. 4

Presentation of the prognostic model in the training set for TCGA-HNSC. (A) Distribution of risk scores for each sample in TCGA-HNSC. (B) Survival status of each sample in the TCGA-HNSC cohort. (C) Heatmap displaying gene expression levels in groups with high- and low-risk scores. (D) Survival curves for groups with high- and low-risk scores. (E) Time-dependent ROC curve. (F) Calibration Curve. TCGA, The Cancer Genome Atlas; HNSC, head and neck squamous cell carcinoma; ROC, receiver operating characteristic

Prognostic model validation

To assess the stability of the model, the same algorithm was used for the GSE41613 dataset, and a score was computed for each sample. Consistent with TCGA cohort findings, patients with high scores exhibited reduced survival times compared to those with lower scores. The validation cohort yielded AUC values of 0.719, 0.666, and 0.691 for 1-, 3-, and 5-year survival, respectively (Fig. 5A-E), with a C-index of 0.64. Furthermore, the calibration curve (Fig. 5F) demonstrated good agreement between predicted and observed survival probabilities, indicating that the model predictions were well-calibrated in the validation cohort.

Additionally, survival analysis of TCGA-HNSC patients was performed using the TNM staging system (AJCC, 7th edition). Kaplan–Meier analysis revealed no statistically significant difference in OS between stages I/II and III/IV (p = 0.307) (Supplementary Fig. S2A). The time-dependent ROC analysis yielded AUC values of 0.672, 0.636, and 0.692 for 1-, 3-, and 5-year survival, respectively, with a C-index of 0.49 (Supplementary Fig. S2B).

Fig. 5

Presentation of the prognostic model in the validation set (GSE41613). (A) Distribution of risk scores for each sample in GSE41613. (B) The survival status of each sample in GSE41613. (C) Heatmap displaying gene expression levels in groups with high- and low-risk scores. (D) Survival curves for groups with high- and low-risk scores. (E) Time-dependent ROC curve. (F) Calibration Curve. TCGA, The Cancer Genome Atlas; HNSC, head and neck squamous cell carcinoma; ROC, receiver operating characteristic

Prognostic significance of the risk score

In the TCGA dataset, univariate Cox regression analysis demonstrated a significant association between the risk score and OS. After adjusting for potential confounders, multivariate analysis confirmed that the risk score is an independent predictor of OS. In the training set, the high-risk group exhibited an HR of 2.55 (95% CI: 1.91–3.39, p = 0.001) compared with the low‐risk group (Fig. 6A). Similarly, in the validation cohort GSE41613, univariate analysis revealed a significant association between the risk score and OS, and multivariate analysis confirmed its independent prognostic value (Fig. 6B). Specifically, the high‐risk group showed an HR of 2.23 (95% CI: 1.25–3.98, p = 0.006) compared with the low‐risk group.

Fig. 6

Independent prognostic value of the score. (A) In the training set TCGA-HNSC, both univariate and multivariate analyses demonstrated that the score is an independent prognostic factor. (B) In the validation set GSE41613, both univariate and multivariate analyses demonstrated that the score is an independent prognostic factor. TCGA, The Cancer Genome Atlas; HNSC, head and neck squamous cell carcinoma

A nomogram was constructed by integrating the score, gender, stage, and smoking status to predict 1-, 3-, and 5-year OS probabilities, with the contribution of each factor to survival risk proportionally represented (Fig. 7A). The calibration curves demonstrated that the nomogram effectively predicted OS rates at 1, 3, and 5 years (Fig. 7B). Decision curve analysis indicated that the nomogram based on the combined model provided superior predictions of patient survival compared to individual prognostic factors (Fig. 7C). As illustrated in Fig. 8, significant differences in risk scores were observed between HPV-positive and HPV-negative groups, as well as among the molecular subtypes of HNSCC patients.

Fig. 7

OS in patients with HNSCC. (A) Nomogram for predicting the probabilities of 1-, 3-, and 5-year OS. (B) Calibration curves suggest that the nomogram demonstrates high accuracy in predicting 1-, 3-, and 5-year OS. (C) DCA shows that the nomogram potentially provides better prediction of OS in patients compared with the use of individual prognostic factors. OS, overall survival; HNSCC, head and neck squamous cell carcinoma; DCA, decision curve analysis

Fig. 8

Comparison of risk scores across HNSCC subgroups defined by clinical features. HNSCC, head and neck squamous cell carcinoma

Examination of immune cell infiltration

The high-risk group exhibited significantly increased infiltration of M0 and M2 macrophages, activated mast cells, resting NK cells, and resting CD4 memory T cells compared to the low-risk group (Fig. 9A). Among the genes in the prognostic model, AREG demonstrated a significant positive correlation with activated mast cells and a negative correlation with resting mast cells. In contrast, ARHGAP4 was positively associated with CD8⁺ T cells, follicular helper T cells, and regulatory T cells (Tregs), while showing a negative correlation with resting memory CD4⁺ T cells (Fig. 9B).

Fig. 9

Analysis of immune cell infiltration. (A) Differences in immune cell infiltration between high- and low-score groups. (B) Correlation between genes in the risk score formula, immune cells, and the score

Analysis of drug sensitivity

Based on drug sensitivity correlation analysis, pazopanib, thapsigargin, imatinib, and NVP-TAE684 were identified as the top four candidate compounds. The log-transformed IC50 values for these agents were significantly lower in the high-risk score group compared to the low-risk group (Fig. 10), indicating that patients in the high-risk group may exhibit increased sensitivity to these drugs.

Fig. 10

Differences in drug sensitivity between high- and low-risk score groups for four drugs

Discussion

HNSCC is a heterogeneous group of malignancies with significant variations in clinical behavior and treatment response. This study examines the role of E2Fs in the molecular heterogeneity of HNSCC. Using consensus clustering of eight E2Fs from the TCGA-HNSC cohort, we identified two distinct subtypes. Notably, despite Cluster 1 showing significantly lower expression of all eight E2F genes, Cluster 2 patients exhibited significantly improved OS. This paradox highlights the complex, context-dependent roles of E2Fs in cancer, where some (e.g., E2F1, E2F3) act as oncogenic drivers while others (e.g., E2F4, E2F5) may serve as tumor suppressors [28,29,30]. This E2F-based subtyping correlates with differences in prognosis, functional enrichments, immune infiltration patterns, and drug sensitivities.

We further identified 211 DEGs between the two clusters. In Cluster 1, upregulated genes were predominantly involved in epidermal development, indicating a more differentiated, squamous-like phenotype, while KEGG analysis highlighted enrichment of the IL-17 signaling pathway, suggesting roles in inflammation and immune responses. IL-17 promotes tumor cell proliferation and survival by activating various signaling pathways, such as the MAPK pathway, thereby affecting patient survival [31]. Elevated IL-17 expression in HNSCC is linked to poor prognosis, as it fosters immune suppression and accelerates tumor progression within the TME [32]. IL-17 promotes tumor progression by disrupting antitumor immune responses through multiple mechanisms, including suppression of chemokine secretion critical for the recruitment of cytotoxic T lymphocytes and NK cells, direct inhibition of their effector functions, and induction of a protumorigenic microenvironment via upregulation of pro-inflammatory cytokines such as IL-6 and TNF-α. Furthermore, IL-17 facilitates immune evasion by promoting M2 macrophage polarization, thereby collectively weakening effective antitumor immunity [33]. Furthermore, IL-17 drives tumor growth and metastasis by activating tumor-associated macrophages and lymphocytes [34]. Abnormal expression of E2F transcription factors is closely associated with cell cycle dysregulation and plays a critical role in tumorigenesis, primarily by accelerating the cell cycle and enhancing DNA replication. Consistent with this role, Cluster 2 showed enrichment in genes related to DNA replication, structural alterations, and double-strand break repair. This suggests that these tumors exhibit higher proliferative activity and altered DNA metabolism. While Cluster 2 was associated with a better overall prognosis, this heightened genomic instability and reliance on DNA repair pathways may increase their sensitivity to DNA-damaging therapies. Overall, these distinct gene expression profiles highlight differences in tumor differentiation, immune response, and proliferation, offering novel insights into the molecular mechanisms of HNSCC and opening avenues for mechanistic studies and targeted therapy development.

Our prognostic model comprises seven genes (AREG, CXCL14, FAM83E, FDCSP, ARHGAP4, EPHX3, and SPINK6), each contributing uniquely to tumor progression and the TME. For example, AREG, a member of the epidermal growth factor family, activates EGFR and its downstream pathways (RAS-RAF-MEK-ERK and PI3K-AKT-mTOR), promoting cancer cell proliferation, migration, and survival. Elevated AREG expression is strongly associated with increased tumor aggressiveness in HNSCC, as it modulates cell cycle–related genes, inhibits apoptosis, and enhances Treg-mediated immune suppression [35,36,37,38]. SPINK6, a serine protease inhibitor, regulates protease activity, especially in skin keratinization. In melanoma, high SPINK6 expression has been linked to tumor aggressiveness and metastasis. EGFR/EphA2 activation promotes tumor cell proliferation and migration [39]. SPINK6 is considered one of the immune regulatory factors associated with Tregs and plays a crucial role in improving the prognosis of HNSCC patients [40]. CXCL14, a highly conserved chemokine primarily expressed in the skin epithelium, plays a multifaceted role in tumor immune evasion by recruiting and maturing immune cells and modulating epithelial motility. Its role appears context-dependent, with higher levels correlating with decreased survival in some cancers [41] and inhibitory effects on tumor proliferation and metastasis in others [42]. Abnormal FAM83E expression has been associated with cancer cell invasion and metastasis, though its mechanisms require further elucidation [43]. FDCSP contributes to immune homeostasis by regulating B-cell maturation and antibody production, thereby influencing tumor immune surveillance and potentially affecting proliferation and metastasis [44]. ARHGAP4, which regulates cytoskeleton reorganization and cell migration, is aberrantly expressed in multiple cancers. Its elevated expression correlates with poor prognosis, as observed in acute myeloid leukemia and colorectal cancer, where it also associates with altered immune cell infiltration [45, 46]. Finally, EPHX3 is involved in metabolic detoxification, and hypomethylation of EPHX3 in patients with oral squamous cell carcinoma is linked to a poor prognosis [47]. Based on correlation analyses between prognostic model genes and immune cells, this study proposes potential therapeutic strategies targeting AREG or ARHGAP4 to modulate the activity of mast cells and T-cell subsets within the TME. For example, inhibiting AREG could effectively reduce mast cell activation, thereby decreasing the subsequent release of pro-tumor factors. Similarly, targeting ARHGAP4 may reduce infiltration or suppressive functions of Tregs, consequently enhancing the activity of anti-tumor CD8⁺ T cells. Although these gene-targeted strategies provide promising therapeutic directions, additional research exploring precise molecular interactions and clinical translation is essential.

The prognostic model developed in this study demonstrated robust predictive performance in both the training cohort and the independent validation dataset, which is comparable to other reported HNSCC prognostic models, with AUC values generally ranging between 0.65 and 0.75 [48]. Unlike some models that rely solely on traditional clinical parameters or a single biomarker [49], our model further integrates data on tumor immune infiltration and drug sensitivity to provide a more comprehensive risk assessment. For TCGA-HNSC, the survival curve analysis based on the TNM staging system did not achieve statistical significance, likely because it fails to fully capture the complexity of patient prognosis. This suggests that the inherent heterogeneity of HNSCC calls for a more robust prognostic model. Although the current LASSO regression effectively selects key genes for risk modeling, it may have limitations in capturing nonlinear features. In the future, more sophisticated machine learning algorithms, such as XGBoost or deep learning approaches, could be considered to enhance the capacity of the model for nonlinear fitting [50].

It is important to note that the risk score from the predictive model may be influenced by platform-specific biases, as the TCGA dataset utilized RNA sequencing while the GSE41613 dataset employed microarray technology. The GSE41613 cohort included 97 samples, whereas the TCGA-HNSC dataset contained 494 samples. Although external validation was reliable, differences in sequencing platforms and sample sizes underscore the need for future platform-independent evaluations to enhance the generalizability and reliability of the prognostic model.

In high-risk HNSCC patients, increased infiltration of immune cells such as M0 and M2 macrophages, activated mast cells, resting NK cells, and resting CD4⁺ memory T cells creates a more immunosuppressive and tumor-promoting TME. Specifically, high levels of M0 macrophages, precursors capable of differentiating into either proinflammatory M1 or anti-inflammatory M2 cells, indicate a reservoir of newly recruited cells poised to respond to TME signals [51]. Upon differentiation, M2 macrophages predominantly support tissue repair, angiogenesis, and matrix remodeling, which facilitate immune evasion and tumor dissemination, key processes linked to a poor prognosis [52]. Elevated activated mast cells further exacerbate tumor progression by releasing mediators, such as histamines and cytokines, that promote angiogenesis and modulate immune responses [53]. These factors collectively shape a TME that is conducive to tumor growth and metastasis. Although NK cells and CD4⁺ memory T cells are typically associated with antitumor activity, their resting state in high-risk patients suggests functional suppression. This suppression is likely driven by inhibitory influences from Tregs or immune checkpoint molecules, such as PD-1/PD-L1, reducing their ability to eliminate tumor cells and contributing to immune tolerance within the TME [54]. This suggests that the immune cells in high-risk patients are functionally suppressed, preventing effective immune surveillance and enhancing tumor survival. It is noteworthy that we analyzed data from the TCGA-HNSC dataset to examine the relationship between HNSCC subtypes constructed based on E2F expression and the immune microenvironment as well as patient outcomes. Our findings indicate that these factors are statistically associated rather than directly causal.

Although tumor purity may influence immune cell estimation, CIBERSORT remains a reliable and widely used method for immune deconvolution. Previous studies in glioblastoma have shown that most immune cell frequencies in tumors are below 0.2%, a range challenging for flow cytometry to accurately measure [55]. Given these complexities, future studies employing single-cell RNA sequencing are crucial for refining cell-type differentiation and elucidating how immune cells interact with tumor cells to affect clinical outcomes. These studies will provide deeper insights into the functional roles of immune cells in the TME and help develop strategies to target immune suppression in high-risk patients.

Our study showed that high-risk patients had lower log(IC50) values for pazopanib, thapsigargin, imatinib, and NVP-TAE684, indicating greater sensitivity to these drugs and suggesting that they may be more effective for treating high-risk patients. Pazopanib, a multi-target tyrosine kinase inhibitor that primarily targets VEGFR, PDGFRs, and c-Kit, has shown variable efficacy in clinical studies both as monotherapy and in combination treatments, particularly in recurrent or metastatic disease [56, 57]. Thapsigargin, a potent inhibitor of the sarcoplasmic/endoplasmic reticulum Ca²⁺ ATPase, disrupts ER calcium homeostasis, leading to ER stress and apoptosis [58]. Although not widely used in HNSCC, its ability to induce apoptosis may benefit cases resistant to conventional therapies. Imatinib, which targets BCR-ABL, c-Kit, and PDGFR, has proven effective in other malignancies such as chronic myeloid leukemia and gastrointestinal stromal tumors; however, its role in HNSCC remains under investigation, with further research needed to identify predictive biomarkers [59, 60]. NVP-TAE684, a small-molecule ALK inhibitor, offers potential as a targeted therapy in HNSCC, although ALK rearrangements are rare in this cancer compared to others like non-small cell lung cancer, necessitating further research to assess its clinical utility [61]. Overall, the observed drug sensitivity in high-risk patients underscores the potential of these agents in personalized HNSCC treatment, warranting further preclinical and clinical validation.

Several limitations must be considered when interpreting our findings. First, the limited HPV status data in the TCGA-HNSC dataset precluded the inclusion of HPV status as a covariate. With only 7 HPV-positive cases and 23 HPV-negative cases, the small sample size hindered robust statistical adjustment. Second, the absence of laboratory experiments limits our mechanistic insights, and the moderate sample size may not fully capture HNSCC heterogeneity, affecting generalizability. Third, the lack of comprehensive clinical validation restricts the applicability of our prognostic model. Finally, our drug sensitivity predictions, derived from GDSC in vitro data and computational modeling, may not directly correlate with clinical outcomes. Both in vitro and in vivo validations are essential to confirm these findings and assess their clinical relevance.

Conclusions

This study provides a comprehensive analysis of the molecular characteristics of HNSCC, emphasizing the pivotal role of E2Fs. By identifying distinct molecular subtypes and constructing a prognostic model, our work enhances the understanding of the molecular landscape in HNSCC. The findings indicate potential biomarkers for diagnosis and novel therapeutic targets that may improve personalized treatment strategies. Future integration of these results with laboratory and clinical research will be crucial for translating them into clinical practice and advancing personalized medicine in HNSCC.