Single-cell analysis of tumor microenvironment immune homeostasis and identification of prognostic biomarkers in head and neck squamous cell carcinoma

Introduction

Head and neck squamous cell carcinoma (HNSCC) is a malignant tumor that originates in areas critical for breathing and swallowing, such as the oral cavity, pharynx, and larynx (1). It is highly invasive and heterogeneous, with cancer-related mortality ranking eighth globally. Owing to the lack of effective early screening tools and molecular biomarkers, 60–70% of patients are already at an advanced stage at diagnosis, and 40–50% present with regional lymph-node metastasis (2,3). The overall 5-year survival rate is approximately 50%, declining to 30–40% in patients with locally advanced disease; recurrence is frequent and the prognosis remains poor (4). At present, treatment for HNSCC is still dominated by surgery, supplemented by conventional modalities such as radiotherapy and chemotherapy; however, these approaches are highly invasive and often fail to ensure patients’ quality of life (5). In the recurrent/metastatic (R/M) setting, immune checkpoint inhibitors (ICIs)—predominantly programmed cell death protein-1 (PD-1)/programmed death-ligand 1 (PD-L1) inhibitors—have become one of the first-line options. In clinical stratification, the combined positive score (CPS) of PD-L1 is commonly used as an important reference, with thresholds (e.g., ≥1 or ≥20) to assess treatment use and potential benefit. Nevertheless, across different clinical scenarios (resectable local disease vs. R/M) and sampling/testing contexts (biopsy vs. whole surgical sections), stratification and efficacy remain variable: different antibody clones and assay platforms have distinct thresholds and concordance; readings from small-volume biopsies and whole surgical sections are often discordant; and in some cohorts the prognostic association of PD-L1 is unstable (6-9). Current evidence indicates that, beyond analytical factors, the spatial and temporal heterogeneity of the tumor microenvironment (TME)—including changes in PD-L1 expression between primary and metastatic sites and before versus after therapy—and its immunosuppressive network constitute key biological sources of instability in PD-L1-based stratification (10,11).

From an etiologic perspective, the development of HNSCC is jointly driven by exogenous risk factors [smoking, heavy alcohol consumption, and high-risk human papillomavirus (HPV)] and intrinsic biological factors (12). As a critical ecological niche for tumorigenesis and progression, the TME exhibits pronounced spatiotemporal heterogeneity and immunosuppressive features. Subtle differences in immune-cell composition, activation states, and cell-cell communication directly influence treatment sensitivity and long-term outcomes (13,14). Within the TME network, T cells are not only the direct targets of PD-1/PD-L1 blockade but also the pivotal link connecting antigen recognition, execution of effector functions, and the formation of immunological memory (15,16). What correlates most closely with treatment outcomes is not a generic “dynamic change”, but rather three measurable states: (I) clonal expansion/contraction (reflecting antigen drive and immunogenicity) (17,18); (II) a differentiation continuum (transition from naïve and memory states to functionally constrained/exhausted states, determining efficacy and durability) (19,20); and (III) intercellular signaling axes (the balance of inhibitory versus costimulatory signals among tumor, myeloid, and stromal cells) (21,22). Mechanism-oriented stratification along these dimensions is more likely to explain inter-individual differences in response and to inform the design of combination strategies. Compared with bulk omics, which “average” information across mixed samples, single-cell transcriptomics not only identifies rare cellular subpopulations and captures functional dynamics that are difficult to observe at the population level, but also systematically resolves intercellular signaling to construct complex communication maps (21,23).

Accordingly, this study integrates transcriptomic sequencing data from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA), using single-cell RNA sequencing (scRNA-seq) data to explore the molecular and dynamic features of the HNSCC TME at the single-cell level. By combining TCGA data with the Cox proportional hazards model, we constructed and validated a prognostic risk prediction model. This model provides new insights for risk stratification and personalized treatment of HNSCC patients, while laying the foundation for the clinical translation of key genes within the tumor immune microenvironment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-999/rc).

Methods

Data download

This study utilized three HNSCC single-cell transcriptomic datasets (GSE172577, GSE180268, and GSE150825) from the GEO database for analysis. The datasets comprised 52 patients with HNSCC and 27 healthy controls, spanning multiple HNSCC subsites (e.g., oral cavity). All single-cell sequencing data were generated using the 10x Genomics platform, which helped reduce batch effects and ensures the comparability of the analysis results. Additionally, RNA sequencing transcriptomic data and corresponding clinical information from 566 HNSCC patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov) for subsequent differential gene analysis and prognostic model construction.

Data quality control

The raw data were stored in FASTQ sequencing file format and aligned to the human reference genome (GRCh38) using the CellRanger software package (v 6.1) to obtain the single-cell gene matrix. In the single-cell gene matrix, cells with fewer than 200 or more than 5,000 RNA molecules, as well as cells with more than 10% mitochondrial gene expression, were considered low-quality cells. The doublets ratio was set to increase by 0.8% for every additional 1,000 cells. Low-quality and double cells were excluded from further analysis (Figures S1,S2).

Clustering analysis

Clustering analysis of the single-cell gene matrix was performed using the R package “Seurat” (v 4.1.1). The “Normalization” and “ScaleData” functions were used for normalization and standardization. Principal components of highly variable genes were computed using the “RunPCA” function. Based on the results of the principal component analysis, important components were identified using the “FindClusters” function with a resolution of 0.5. Uniform manifold approximation and projection (UMAP) was performed using RunUMAP, and cell clusters were projected into two-dimensional space for visualization with UMAPPlot (24). Finally, the “FindAllMarkers” function was used to identify marker genes within the clusters.

Cell-cell interaction analysis

Cell-cell interaction analysis was performed using the R package “CellChat” (v5.0.1). The communication probabilities between cells in the analysis were computed to examine the interactions within the dataset. Subsequently, the “netVisual_circle” function was used to visualize the cell interactions based on the number and weight of the interactions. Ligand-receptor interactions between cell types were predicted based on the expression of ligands and receptors. During the analysis, only ligands and receptors expressed in more than 10% of the cells within the clusters were considered.

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

GO and KEGG enrichment analyses were performed using the “enrichGO” and “enrichKEGG” functions from the R package “clusterProfiler” (v 4.6.2) (25). Results were considered statistically significant when the P value and false discovery rate (FDR) were both below 0.05. Genes that met these criteria were used for GO and KEGG enrichment analysis. Significant enrichment terms related to specific biological functions, processes, or metabolic pathways were identified from the gene sets.

Pseudotime analysis

Pseudotime analysis was performed using the R package monocle3 (v1.4.25). The matrix used in Monocle3 was derived by converting the “Seurat” object into a cell_data_set object. Genes used for ordering were selected from differentially expressed genes (DEGs) between clusters defined in the clustering analysis. The “reduce_dimension” and “order_cells” functions were used to order the cells, and the “plot_cells” function was employed to visualize the cell trajectories, selecting an initial cell state and pseudotime direction. The “graph_test” function was used to identify DEGs along the cell trajectory, and the top ten genes were selected for visualization based on their Moran’s I values.

Identification of DEGs in TCGA-HNSC

The “DESeq2” (v 1.48.1) package was used to analyze the transcriptomic data from the case and normal control groups in the TCGA database, identifying DEGs. DEGs with a log₂ fold change >2 and P<0.05 were defined as upregulated genes, while DEGs with a log₂ fold change <−2 and P<0.05 were defined as downregulated genes. Genes with values between these thresholds were classified as stable genes.

Prognostic prediction model construction

To identify genes significantly associated with the prognosis of HNSCC patients, a univariate Cox regression model was first applied to preliminarily screen all candidate genes, with a P value threshold of <0.05 to determine significantly associated genes. Next, a least absolute shrinkage and selection operator (LASSO) regression model was used for further selection and compression of these genes, with the optimal penalty parameters determined through 10-fold cross-validation (lambda.min =0.0184 and lambda.1se =0.1080). Genes with non-zero regression coefficients from the LASSO regression were used to construct a multivariate Cox regression model. In the multivariate Cox regression model, the independent effects of each gene on overall patient survival were assessed, and regression coefficients, hazard ratios (HRs), 95% confidence intervals (CIs), and P values were extracted to determine their prognostic value. Based on the calculated risk scores from the model, patients were stratified into high-risk and low-risk groups according to the median risk score. The discriminatory ability of these groups in survival prediction was further validated.

Statistical analysis

All analyses were conducted in R (v4.4.1). Unless otherwise specified, tests were two-sided, multiple comparisons were controlled using the Benjamini-Hochberg FDR, and statistical significance was defined as P<0.05 or FDR <0.05. P values are reported in accordance with the journal’s formatting policy. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Results

Identification of tumor-infiltrating cell types in human HNSCC based on 10x scRNA-seq

After removing doublets and performing quality control, a total of 202,411 cells from the single-cell datasets were retained for subsequent analysis. Among these, 145,828 cells were derived from 52 HNSCC patients, and 56,583 cells were from 27 healthy controls. Following batch effect correction, each cell was mapped to a comparable space in the UMAP projection, resulting in the identification of 22 clusters. Based on marker genes curated from published literature and databases (26,27), clusters were manually annotated, yielding 11 cell types, including stem cells, T cells, B cells, myeloid cells, and tumor cells (Figure 1A-1C). The distribution of cell types in the clusters revealed that stem cells comprised the largest proportion of all cell types, followed by T cells, B cells, and myeloid cells (Figure 1D,1E).

Figure 1 Identification of cell types in HNSCC single-cell transcriptomic data. (A) Clustering of HNSCC cell clusters; (B) HNSCC cell type identification plot; (C) HNSCC gene expression dot plot; (D) comparison of cell types between HNSCC patient and control groups; (E) proportions of different cell types in HNSCC. HNSCC, head and neck squamous cell carcinoma; NK, natural killer; TD, T cell differentiation state; UMAP, uniform manifold approximation and projection.

Elucidation of intercellular communication in HNSCC through ligand-receptor network analysis

After establishing the cell-type framework and relative composition, we further evaluated ligand-receptor interactions among distinct cellular populations to delineate upstream regulatory networks that may shape immune states. This study utilized the CellChat package to identify ligand-receptor interactions and cell communication between major cell types. The analysis was conducted using the CellChatDB.human database, which includes 61.8% paracrine/autocrine signaling interactions, 21.7% extracellular matrix (ECM) receptor interactions, and 16.8% cell-cell communication interactions (Figure 2A). The relevant cell-cell interaction network showed that T cells exhibited high degree/strength, with dense interactions with tumor cells, macrophages, and stromal cells; stem cells were also tightly connected with macrophages, T cells, and epithelial cells, suggesting possible involvement in immunosuppressive or repair-related regulation within the microenvironment (28) (Figure 2B,2C). At the molecular level, we identified multiple significant pairs related to human leukocyte antigen (HLA) molecules within the communication network (including HLA-E/HLA-F-CD8A/CD8B), suggesting a potential nonclassical major histocompatibility complex (MHC) class I-centered communication axis between tumor-associated cell populations and CD8⁺ T cells (Figure 2D). HLA molecules are key components of the MHC and are commonly associated with immune regulation (29). High expression of HLA-E and HLA-F may, through their interactions with CD8⁺ T cells, induce T-cell exhaustion and thereby suppress antitumor immune responses within the tumor (30,31).

Figure 2 Cell interaction network. (A) Diagram of the CellChatDB.human database; (B) overall diagram of interactions between different cell types; (C) overall diagram of interaction weights between different cell types; (D) heatmap of receptor-ligand pair expression for cell type interactions. ECM, extracellular matrix; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Mapping T cell clusters and their developmental trajectories in HNSCC

Given the centrality of T cells in the communication network, we further explored the functional characteristics of T cells in the HNSCC TME and their role in immune evasion. Based on scRNA-seq data from patients with HNSCC, a total of 21,671 T cells were extracted. After UMAP reduction and clustering, T cells were partitioned into the major subsets of exhausted T cells (T_ex), memory T cells (T_mem), regulatory T cells (T_reg), and naive T cells (T_n) (Figure 3A,3B). In terms of composition, T_mem constituted the largest proportion, followed by T_reg (Figure 3C).

Figure 3 Clustering analysis of T cell subpopulations in HNSCC. (A) Clustering analysis of T cell subpopulations in the patient group; (B) T cell subtype identification; (C) proportion of T cell subtypes in each cell cluster; enrichment analysis: (D) GO enrichment analysis for BPs, CCs, and MFs; (E) KEGG enrichment analysis; (F) Pseudotime analysis of T cell subtypes—cell ordering plot; (G) pseudotime gene expression plot for T cell subtypes. BP, biological process; CC, cellular component; GO, Gene Ontology; HNSCC, head and neck squamous cell carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; UMAP, uniform manifold approximation and projection.

Pathway enrichment analysis of T-cell subsets showed that, in KEGG terms, enrichment was prominent in the T-cell receptor (TCR) signaling pathway, antigen processing and presentation, T helper 17 (Th17) cell differentiation, leukocyte transendothelial migration, and Proteoglycans in cancer. For GO, biological processes were enriched for T-cell differentiation, cell-cell adhesion, and antiviral-related processes; molecular functions were enriched for transcription factor binding and small GTPase-related activity; and cellular components were enriched for nuclear-associated structures/processes (e.g., RNA splicing and chromatin-related) (Figure 3D,3E). Monocle3 was used for pseudotime analysis to depict the continuum of states. The trajectory progressed from early-pseudotime Tn to Tmem and subsequently bifurcated at late pseudotime toward Tex and Treg (Figure 3F). The top ten genes most strongly associated with the trajectory, ranked by Moran’s I, displayed marked dynamics along pseudotime: CCL5 and FOXP3 increased at late pseudotime; effector-related genes GZMK and NKG7 were higher at early pseudotime and then declined; metabolism/signaling-related genes MTRNR2L12 and FXYD6 varied in a time-dependent manner; and TNFRSF4 (OX40) showed trajectory-specific expression patterns (Figure 3G).

Key molecules identified from TCGA-HNSC data and their clinical significance

The state continuity at the single-cell level and the candidate regulatory genes suggested potential functional nodes. To test their robustness and clinical relevance at the population level, we performed differential expression analysis in TCGA-HNSC and intersected the results with single-cell marker genes to obtain a candidate set consistent across data types. Differential expression analysis of TCGA-HNSC data identified a total of 3,656 DEGs, including 2,052 upregulated genes and 1,604 downregulated genes (tumor samples n=522, control samples n=44). GO enrichment analysis revealed that the upregulated genes were primarily involved in receptor-ligand activity, DNA-binding transcription factor activity, and RNA polymerase II-specific transcriptional regulation, reflecting enhanced transcriptional activation and intercellular communication functions within the TME. In contrast, downregulated genes were enriched in actin binding, passive transmembrane transport activity, and channel activity, indicating impaired cytoskeletal dynamics and metabolic regulation in tumor tissues. These results are consistent with the known molecular dysfunctions in HNSCC, supporting the reliability of the findings (27,32) (Figure 4A,4B).

Figure 4 Overview and identification of key risk genes in TCGA HNSCC data analysis. (A) Volcano plot of DEGs in TCGA-HNSC dataset; (B) GO enrichment analysis of upregulated genes in the DEGs; (C) Venn diagram showing overlap between upregulated genes in DEGs and marker genes from single-cell data; (D) heatmap of survival contributions of overlapping genes in different tumor types; (E) survival curve for low-risk genes in HNSCC; (F) survival curve for high-risk genes in HNSCC. DEGs, differentially expressed genes; GO, Gene Ontology; HNSCC, head and neck squamous cell carcinoma; HR, hazard ratio; TCGA, The Cancer Genome Atlas.

After intersecting the upregulated DEGs with the marker genes identified from the single-cell data, 65 overlapping genes were obtained (Figure 4C). These overlapping genes were further analyzed for their prognostic value using the gene expression profiling interactive analysis 2 (GEPIA2) tool, focusing on 10 genes that play a critical regulatory role in T cell dynamic development. Survival analysis identified 15 risk genes in total, among which eight high-risk genes (e.g., TGFBI, PDPN, PLAU) with high expression were significantly associated with poorer overall survival (OS), whereas seven low-risk genes (e.g., FOXP3, TNFRSF4, NKG7) with high expression were associated with better OS (Figure 4D). Kaplan-Meier curve analysis indicated that patients with a high-risk gene profile had significantly lower survival rates compared to those with a low-risk gene profile (log-rank P<0.001), emphasizing the importance of these genes as potential prognostic biomarkers (Figure 4E,4F). These high-risk genes may act as oncogenic drivers that promote tumor progression and microenvironment remodeling, while the low-risk genes may exert tumor-suppressive effects by enhancing immune responses. These findings provide a solid foundation for further exploring the functional roles of these genes in HNSCC and their potential as therapeutic targets or prognostic indicators.

Construction and validation of the gene prognostic model

After obtaining candidate genes that were upregulated in tumor tissues and consistent with single-cell features, a six-gene core set (SERPINH1, PLAU, INHBA, TNFRSF4, CXCL13, STAG3) was selected from seven low-risk and eight high-risk genes using univariate Cox regression and LASSO regression. A prognostic risk model was constructed based on a multivariate Cox regression model, and the corresponding risk scores were calculated. The TCGA-HNSC tumor samples were divided into high-risk and low-risk groups based on the median risk score. Among the six genes, SERPINH1, PLAU, and INHBA were significantly expressed in the high-risk group, while TNFRSF4, CXCL13, and STAG3 were significantly expressed in the low-risk group (Figure 5A).

Figure 5 Screening of risk genes for prognostic prediction in HNSCC patients. (A) Box plot of expression levels of model-related genes between the two groups; (B) forest plot of clinical factors related to OS identified by univariate Cox regression analysis; (C) Kaplan-Meier survival curve based on prognostic risk score; (D) dynamic prediction performance ROC curve for OS at 1-, 3-, and 5-year time points. Asterisks denote statistical significance: ****, P<0.0001. AUC, area under the curve; CI, confidence interval; HNSCC, head and neck squamous cell carcinoma; OS, overall survival; ROC, receiver operating characteristic.

Subsequently, we used univariate Cox regression analysis in the TCGA cohort to investigate whether the risk score independently predicts patient prognosis, regardless of traditional clinical factors. The results showed that the risk score (HR: 2.40; 95% CI: 1.7–3.4; P<0.001) is an independent predictor of OS in HNSCC patients (Figure 5B). Furthermore, Kaplan-Meier survival curves and log-rank tests revealed a significant difference in OS between high-risk and low-risk HNSCC patients (Figure 5C). The risk score demonstrated particularly strong predictive performance in the TCGA cohort, with areas under the curve (AUCs) of 0.640, 0.668, and 0.596 for 1-, 3-, and 5-year OS predictions, respectively, indicating its prognostic potential (Figure 5D).

Discussion

Tumor cells’ ability to adapt to or evade immune surveillance is one of the hallmark features of cancer, which leads to significant variability in patient responses to ICIs (33). With the continuous advancement in TME research, anti-tumor immunotherapy, which activates the body’s immune system, has emerged as the fourth major treatment modality following surgery, radiotherapy, and chemotherapy (34). T cells, as the central players in tumor immune responses, play a critical role in tumor initiation and progression. Their differentiation and functional status within the TME directly influence the effectiveness of immune therapy (35). In this study, based on scRNA-seq data from HNSCC, we systematically analyzed the high heterogeneity of T cells within the TME, identifying multiple T cell subtypes, including T_ex and T_reg. We further utilized pseudotime analysis to reveal the dynamic differentiation trajectories of T cells within the TME, from their initial phenotype to exhausted and immune-suppressive phenotypes. This process is regulated by key genes such as CCL5, FOXP3, and NKG7. Notably, high expression of CCL5 is closely associated with the dysfunction of exhausted T cells and has been reported to contribute to immune suppression in various solid tumors by promoting the aggregation and functional impairment of exhausted T cells (36,37). The expression of FOXP3, a significant marker of regulatory T cell immune-suppressive function, suggests that these genes play a vital role in tumor immune evasion (38).

Previous studies have highlighted the importance of genetic characteristics in guiding cancer treatment and prognostic evaluation (39,40). In this study, we further intersected the single-cell marker genes with the DEGs from TCGA-HNSC data, identifying 15 genes significantly associated with prognosis. Based on LASSO and multivariate Cox regression models, we selected six key genes (SERPINH1, PLAU, INHBA, TNFRSF4, CXCL13, and STAG3) with significant prognostic potential and constructed an accurate prognostic risk model. Integrating scRNA-seq with phenotypic data from the TCGA database, we identified gene sets with opposite directions with respect to prognosis: the SERPINH1-PLAU-INHBA set was associated with poorer outcomes, whereas the TNFRSF4-CXCL13-STAG3 set was associated with better outcomes. Among the high-risk genes, SERPINH1 (HSP47) is related to collagen maturation, ECM densification, and epithelial-mesenchymal transition (EMT), commonly observed in reactive stroma and concomitant with EMT phenotypes (41); PLAU [encoding urokinase-type plasminogen activator (uPA)] mediates fibrinolysis and degradation of the basement membrane and ECM, and the uPA–uPAR axis activates the PI3K/AKT and MAPK signaling pathways (42); INHBA is linked to the activin-TGF-β/cancer-associated fibroblast (CAF) axis, consistent with reports of invasive and immune-excluded phenotypes (43). In contrast, among the low-risk genes, TNFRSF4 (OX40) points to T-cell costimulation and memory maintenance and can attenuate Treg suppression in specific contexts (44); CXCL13 is associated with B-cell recruitment and intratumoral formation of tertiary lymphoid structures (TLS), often seen in an immune-activated microenvironment (45); STAG3, a component of the cohesin complex, participates in maintaining genome stability, consistent with reports of better outcomes (46). Based on the concordant associations across single-cell and TCGA analyses, we infer that the low-risk group tends toward an immune-activated microenvironment, whereas the high-risk group tends toward an ECM/fibrinolysis- and TGF-β-dominant microenvironment. It is noteworthy that FOXP3 was associated with better survival in univariate analysis but was no longer independent in the multigene model; this signal may reflect co-occurrence and collinearity with immune-activation/TLS markers (e.g., CXCL13) rather than a direct effect of Treg-dominated suppression (47,48).

From a clinical perspective, the six-gene risk model in this study showed moderate discrimination for 3-year overall survival (AUC =0.668). In the absence of integration of key clinical variables (e.g., HPV/p16 status and treatment modalities), it can be regarded as a gene-based baseline risk indicator. Although treatment heterogeneity and follow-up bias in TCGA, as well as limited clinical annotation in public datasets, may constrain model performance, the score still showed independent prognostic value and significant survival separation in this cohort, suggesting that it can serve as a complementary dimension to clinical information. Future work may develop a nomogram in an independent cohort with complete clinical data by combining age, stage, HPV/p16, and treatment variables, and systematically evaluate clinical utility using time-dependent AUC, concordance index (C-index), calibration curves, and decision curve analysis (DCA). On this basis, the model could be used for postoperative follow-up stratification, optimization of therapeutic strategies, and trial enrollment enrichment, thereby providing practical value in reducing overtreatment and improving resource allocation.

Strengths and limitations

However, despite the systematic analysis of the HNSCC microenvironment and prognostic features through the integration of scRNA-seq and bulk RNA-seq data from TCGA, there are still certain limitations in this study: (I) the GEO and TCGA datasets may not fully represent the diverse HNSCC patient population, and sample selection bias may affect the generalizability of the results. The proposed prognostic model and biomarkers need validation in larger cohorts. (II) Although key genes were identified, their specific biological functions have not yet been validated through in vitro or in vivo experiments. Future studies can combine gene knockdown/overexpression in vitro and in vivo animal models to further investigate the roles of genes such as SERPINH1, PLAU, and INHBA in HNSCC invasion, metastasis, and immune regulation. Additionally, the interactions between TNFRSF4, CXCL13, STAG3, and other immune checkpoints (e.g., CTLA-4, LAG-3, TIM-3) need further exploration to assess their clinical feasibility and safety in personalized therapy. By integrating the expression profiles of these six core genes with patient factors such as HPV infection status, staging, and treatment regimens, a more precise prognostic assessment system can be developed, providing stronger molecular support for individualized diagnostic and therapeutic strategies for HNSCC patients.

Conclusions

In summary, this study systematically revealed the heterogeneity and dynamic changes of T cells within the TME by integrating single-cell data with TCGA-HNSC data. A prognostic risk model based on key genes was constructed, enabling precise risk stratification for HNSCC patients. The findings provide a theoretical basis for immune therapy and precision prognostic assessment in HNSCC and lay the foundation for gene-based personalized treatment strategies. In the future, further validation through additional independent cohorts and in-depth functional experiments will likely advance the development of precision treatment and prognostic management for HNSCC.

Acknowledgments

None.

Footnote

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

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

Funding: This work was supported by the National Natural Science Foundation of China (Nos. 32071462 and 32471520).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-999/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. The 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. Johnson DE, Burtness B, Leemans CR, et al. Head and neck squamous cell carcinoma. Nat Rev Dis Primers 2020;6:92. [Crossref] [PubMed]
2. Singh D, Vignat J, Lorenzoni V, et al. Global estimates of incidence and mortality of cervical cancer in 2020: a baseline analysis of the WHO Global Cervical Cancer Elimination Initiative. Lancet Glob Health 2023;11:e197-206. [Crossref] [PubMed]
3. von Witzleben A, Wang C, Laban S, et al. HNSCC: Tumour Antigens and Their Targeting by Immunotherapy. Cells 2020;9:2103. [Crossref] [PubMed]
4. Kawazoe A, Fukuoka S, Nakamura Y, et al. Lenvatinib plus pembrolizumab in patients with advanced gastric cancer in the first-line or second-line setting (EPOC1706): an open-label, single-arm, phase 2 trial. Lancet Oncol 2020;21:1057-65. [Crossref] [PubMed]
5. Machiels JP, René Leemans C, Golusinski W, et al. Squamous cell carcinoma of the oral cavity, larynx, oropharynx and hypopharynx: EHNS-ESMO-ESTRO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Ann Oncol 2020;31:1462-75. [Crossref] [PubMed]
6. Girolami I, Pantanowitz L, Barberis M, et al. Challenges facing pathologists evaluating PD-L1 in head & neck squamous cell carcinoma. J Oral Pathol Med 2021;50:864-73. [Crossref] [PubMed]
7. Nocini R, Vianini M, Girolami I, et al. PD-L1 in oral squamous cell carcinoma: A key biomarker from the laboratory to the bedside. Clin Exp Dent Res 2022;8:690-8. [Crossref] [PubMed]
8. Cerbelli B, Girolami I, Eccher A, et al. Evaluating programmed death-ligand 1 (PD-L1) in head and neck squamous cell carcinoma: concordance between the 22C3 PharmDx assay and the SP263 assay on whole sections from a multicentre study. Histopathology 2022;80:397-406. [Crossref] [PubMed]
9. Wang C, Hahn E, Slodkowska E, et al. Reproducibility of PD-L1 immunohistochemistry interpretation across various types of genitourinary and head/neck carcinomas, antibody clones, and tissue types. Hum Pathol 2018;82:131-9. [Crossref] [PubMed]
10. Bill R, Faquin WC, Pai SI. Assessing PD-L1 Expression in Head and Neck Squamous Cell Carcinoma: Trials and Tribulations. Head Neck Pathol 2023;17:969-75. [Crossref] [PubMed]
11. Lin X, Kang K, Chen P, et al. Regulatory mechanisms of PD-1/PD-L1 in cancers. Mol Cancer 2024;23:108. [Crossref] [PubMed]
12. Duan Q, Zhang H, Zheng J, et al. Turning Cold into Hot: Firing up the Tumor Microenvironment. Trends Cancer 2020;6:605-18. [Crossref] [PubMed]
13. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature 2017;541:321-30. [Crossref] [PubMed]
14. Bruni D, Angell HK, Galon J. The immune contexture and Immunoscore in cancer prognosis and therapeutic efficacy. Nat Rev Cancer 2020;20:662-80. [Crossref] [PubMed]
15. Wei SC, Duffy CR, Allison JP. Fundamental Mechanisms of Immune Checkpoint Blockade Therapy. Cancer Discov 2018;8:1069-86. [Crossref] [PubMed]
16. Pardoll DM. The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer 2012;12:252-64. [Crossref] [PubMed]
17. Wang D, Fang J, Wen S, et al. A comprehensive profile of TCF1+ progenitor and TCF1- terminally exhausted PD-1+CD8+ T cells in head and neck squamous cell carcinoma: implications for prognosis and immunotherapy. Int J Oral Sci 2022;14:8. [Crossref] [PubMed]
18. La Manno G, Soldatov R, Zeisel A, et al. RNA velocity of single cells. Nature 2018;560:494-8. [Crossref] [PubMed]
19. Im SJ, Hashimoto M, Gerner MY, et al. Defining CD8+ T cells that provide the proliferative burst after PD-1 therapy. Nature 2016;537:417-21. [Crossref] [PubMed]
20. Miller BC, Sen DR, Al Abosy R, et al. Subsets of exhausted CD8(+) T cells differentially mediate tumor control and respond to checkpoint blockade. Nat Immunol 2019;20:326-36. [Crossref] [PubMed]
21. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity. Nat Rev Immunol 2018;18:35-45. [Crossref] [PubMed]
22. Binnewies M, Roberts EW, Kersten K, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med 2018;24:541-50. [Crossref] [PubMed]
23. Riaz N, Havel JJ, Makarov V, et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 2017;171:934-949.e16. [Crossref] [PubMed]
24. Lin P, Troup M, Ho JW. CIDR: Ultrafast and accurate clustering through imputation for single-cell RNA-seq data. Genome Biol 2017;18:59. [Crossref] [PubMed]
25. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
26. Hu C, Li T, Xu Y, et al. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res 2023;51:D870-6. [Crossref] [PubMed]
27. Puram SV, Tirosh I, Parikh AS, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 2017;171:1611-1624.e24. [Crossref] [PubMed]
28. 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]
29. Liu B, Shao Y, Fu R. Current research status of HLA in immune-related diseases. Immun Inflamm Dis 2021;9:340-50. [Crossref] [PubMed]
30. Li Y, Li Z, Tang Y, et al. Unlocking the therapeutic potential of the NKG2A-HLA-E immune checkpoint pathway in T cells and NK cells for cancer immunotherapy. J Immunother Cancer 2024;12:e009934. [Crossref] [PubMed]
31. Garcia-Beltran WF, Hölzemer A, Martrus G, et al. Open conformers of HLA-F are high-affinity ligands of the activating NK-cell receptor KIR3DS1. Nat Immunol 2016;17:1067-74. [Crossref] [PubMed]
32. Datta A, Deng S, Gopal V, et al. Cytoskeletal Dynamics in Epithelial-Mesenchymal Transition: Insights into Therapeutic Targets for Cancer Metastasis. Cancers (Basel) 2021;13:1882. [Crossref] [PubMed]
33. Agudo J, Park ES, Rose SA, et al. Quiescent Tissue Stem Cells Evade Immune Surveillance. Immunity 2018;48:271-285.e5. [Crossref] [PubMed]
34. Cai XJ, Zhang HY, Zhang JY, et al. Bibliometric analysis of immunotherapy for head and neck squamous cell carcinoma. J Dent Sci 2023;18:872-82. [Crossref] [PubMed]
35. Heft Neal ME, Spielbauer KK, Spector ME. Time to Surgery and Survival in Head and Neck Cancer. Ann Surg Oncol 2021;28:602-3. [Crossref] [PubMed]
36. Li C, Chen S, Liu C, et al. CCR5 as a prognostic biomarker correlated with immune infiltrates in head and neck squamous cell carcinoma by bioinformatic study. Hereditas 2022;159:37. [Crossref] [PubMed]
37. Aldinucci D, Borghese C, Casagrande N. The CCL5/CCR5 Axis in Cancer Progression. Cancers (Basel) 2020;12:1765. [Crossref] [PubMed]
38. Sakaguchi S, Miyara M, Costantino CM, et al. FOXP3+ regulatory T cells in the human immune system. Nat Rev Immunol 2010;10:490-500. [Crossref] [PubMed]
39. Yang C, Cheng X, Gao S, et al. Integrating bulk and single-cell data to predict the prognosis and identify the immune landscape in HNSCC. J Cell Mol Med 2024;28:e18009. [Crossref] [PubMed]
40. Wang SY, Dang W, Richman I, et al. Cost-Effectiveness Analyses of the 21-Gene Assay in Breast Cancer: Systematic Review and Critical Appraisal. J Clin Oncol 2018;36:1619-27. [Crossref] [PubMed]
41. Murugesan D, Kannan B, As SG, et al. Alteration of SERPINH1 is associated with elevated expression in head and neck squamous cell carcinomas and its clinicopathological significance. J Stomatol Oral Maxillofac Surg 2024;125:101811. [Crossref] [PubMed]
42. Lallemant B, Evrard A, Combescure C, et al. Clinical relevance of nine transcriptional molecular markers for the diagnosis of head and neck squamous cell carcinoma in tissue and saliva rinse. BMC Cancer 2009;9:370. [Crossref] [PubMed]
43. Shimizu S, Seki N, Sugimoto T, et al. Identification of molecular targets in head and neck squamous cell carcinomas based on genome-wide gene expression profiling. Oncol Rep 2007;18:1489-97.
44. Vetto JT, Lum S, Morris A, et al. Presence of the T-cell activation marker OX-40 on tumor infiltrating lymphocytes and draining lymph node cells from patients with melanoma and head and neck cancers. Am J Surg 1997;174:258-65. [Crossref] [PubMed]
45. Vachani A, Nebozhyn M, Singhal S, et al. A 10-gene classifier for distinguishing head and neck squamous cell carcinoma and lung squamous cell carcinoma. Clin Cancer Res 2007;13:2905-15. [Crossref] [PubMed]
46. Tian G, Fu Y, Zhang D, et al. Identification of four key prognostic genes and three potential drugs in human papillomavirus negative head and neck squamous cell carcinoma. Cancer Cell Int 2021;21:167. [Crossref] [PubMed]
47. Shang B, Liu Y, Jiang SJ, et al. Prognostic value of tumor-infiltrating FoxP3+ regulatory T cells in cancers: a systematic review and meta-analysis. Sci Rep 2015;5:15179. [Crossref] [PubMed]
48. Cabrita R, Lauss M, Sanna A, et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. Nature 2020;577:561-5. [Crossref] [PubMed]