The TCGA-HNSC cohort, encompassing 502 HNSC and 44 normal samples, was retrieved from the TCGA database (https://portal.gdc.cancer.gov/), with 501 HNSC samples possessing complete clinical data included in subsequent analyses. These samples were randomized into training and test risk groups (7:3 ratio) using the cart R package (https://cran.r-project.org/web/packages/caret/index.html), drawing on relevant clinical information. Additionally, 40 patients with comprehensive follow-up data were sourced from the ICGC as an external validation set.

Differential expression analysis between normal and tumor groups was meticulously executed using the ‘limma’ R package [26, 27, 28], adhering to stringent criteria of an absolute log fold change ($|$logFC$|$) greater than 1 and an adjusted p-value threshold below 0.05. This approach facilitated the identification of genes significantly altered in HNSCC compared to normal tissue. Further, we employed the ‘clusterProfiler’ R package (version 4.1.3) [29, 30] to conduct comprehensive enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, setting a significance threshold at a p-value of less than 0.05. The GO enrichment analysis delved into the molecular functions (MF), biological processes (BP), and cellular components (CC) associated with these differentially expressed genes. Concurrently, the KEGG pathway analysis provided insights into the biological pathways implicated in HNSCC, thereby enhancing our understanding of the molecular underpinnings of this complex disease.

TCGA-HNSC cohorts were segregated into distinct groups based on the consensus expression of soluble media-associated differentially expressed genes (DE-SMRGs) using ‘ConsensionClusterPlus’ R software (Version 4.1.3, Chapel Hill, NC, USA) [31]. The MCPcounter R package (Version 4.1.3, Paris, France) [32] was employed to evaluate the content of immune-related cells in these clusters and to generate violin plots.

The prognostic value of SMRGs in the HNSCC cohort was meticulously evaluated through univariate Cox regression analysis. This statistical approach allowed us to determine the individual impact of each SMRG on patient survival. Further, to refine our analysis and identify the most influential genes within the vast array of SMRGs, we employed LASSO Cox regression analysis using the “glmnet” R package [33]. This method effectively pinpointed key genes and their corresponding regression coefficients, reducing the risk of overfitting by penalizing the complexity of the model. The risk scores for each patient were then calculated using a comprehensive formula that combines the expression levels of these key mRNAs and their respective regression coefficients. The formula is structured as follows: risk score = Expression_mRNA1 $\times$ Coef_mRNA1 + Expression_mRNA2 $\times$ Coef_mRNA2 + … + Expression_mRNAn $\times$ Coef_mRNAn. This risk scoring system encapsulates the cumulative impact of multiple gene expressions on the prognosis of HNSCC, enabling a more nuanced and precise prediction of clinical outcomes.

Risk scores for patients in our study were accurately computed using the “survminer” R package, based on the established model equations derived from the combination of gene expressions and their corresponding coefficients. This calculation enabled the categorization of patients into high-risk (HR) and low-risk (LR) groups, a division based on the median value of the calculated risk scores. Following this stratification, survival curves for each group were meticulously plotted to visually represent and compare the survival probabilities over time between these two distinct risk categories. To quantitatively assess the predictive accuracy of our prognostic model, the concordance index (C-index) was determined using the ‘pec’ R program. The C-index is a crucial metric in survival analysis, providing a measure of the model’s ability to accurately predict survival outcomes [34]. Additionally, we utilized the ‘survivalROC’ R package for conducting time-dependent Receiver Operating Characteristic (ROC) curve analysis. This analysis is essential for evaluating the predictive potential of the genetic markers identified in our study over time, offering insights into the temporal robustness and reliability of the prognostic model in assessing patient outcomes in HNSCC.

The evaluation of the risk score’s effectiveness as an independent prognostic factor for HNSCC was thoroughly conducted using both univariate and multivariate Cox regression models. These models are instrumental in discerning the extent to which the calculated risk score can predict patient outcomes independently of other clinical variables. Furthermore, to translate our findings into a practical tool for clinical use, we utilized the “rms” R package [35] to create detailed nomograms. These nomograms integrate several key parameters, including the risk score derived from our model, patient age, tumor stage, and the expression levels of the model genes. They are designed to predict the 1, 3, and 5-year overall survival rates for patients in the TCGA HNSC cohort.

We employed a multifaceted approach to quantify immune infiltration scores, utilizing a suite of established methodologies. These included XCELL [36, 37], TIMER [38, 39], QUANTISEQ [40, 41], MCPCOUNT [32], EPIC [42], CIBERSORT [43, 44]and CIBERSORT-ABS [45]. The integration of these diverse tools enabled a thorough analysis of the immune landscape within the HNSCC tumor microenvironment. We conducted Spearman correlation analysis to investigate the relationship between the calculated risk scores and the presence of various immune cells, providing insights into the interplay between genetic risk factors and the immune milieu. Furthermore, the single-sample Gene Set Enrichment Analysis (ssGSEA) method was applied to determine the immune cell characteristics in HNSCC patients. This approach was instrumental in differentiating the immune infiltration status between high-risk (HR) and low-risk (LR) groups, offering a nuanced understanding of how risk levels correlate with immune response.

Our analysis extended to examining 20 suppressive immune checkpoints identified in Auslander’s study [46], focusing on their variations between the HR and LR groups. Additionally, we referenced a gene set pertaining to the cancer immune cycle from Xu et al.’s website [44] (http://biocc.hrbmu.edu.cn/TIP/) and drew upon Mariathasan’s research for a list of genes associated with favorable responses to the anti-PD-L1 drug atezolizumab [47]. To assess the enrichment of gene features related to the cancer immune cycle and immunotherapy response, we employed the Gene Set Variation Analysis (GSVA) method [48] was employed. A p-value of less than 0.05 was considered indicative of a significant difference between the HR and LR groups. The ‘ggcor’ R package (version 4.1.3) was utilized to further analyze the association between risk scores and these genetic characteristics, providing a comprehensive overview of how genetic risk factors might influence the efficacy of immunotherapeutic interventions in HNSCC.

The evaluation of treatment response in high-risk (HR) and low-risk (LR) patient groups with HNSCC was meticulously conducted using half-maximal inhibitory concentration (IC50) data. This data was sourced from the Genomics of Drug Sensitivity in Cancer (GDSC) database, accessible at https://www.cancerrxgene.org/. Utilizing the “pRRophetic” R package (version 4.1.3) [49], we were able to predict and compare the response of HR and LR groups to various anticancer drugs.

The Tumor Immune Single-Cell Hub (TISCH) database, accessible at http://tisch.comp-genomics.org, played a pivotal role in our research, particularly in the study of tumor microenvironment (TME) heterogeneity in HNSCC. As a comprehensive repository of single-cell RNA sequencing data, TISCH specializes in providing detailed insights into the TME across a variety of cancer types, datasets, and cell types [50].

The verification of the expression levels of soluble mediator-related genes (SMGs) in our study on HNSCC was significantly bolstered by utilizing immunohistochemical sections from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). HPA is an extensive resource that integrates proteomic, transcriptomic, and systems biology data to provide a comprehensive view of human protein expression across various tissues and cells.

In our comprehensive study, all statistical analyses were rigorously performed using R software, version 4.1.3. We adhered to stringent statistical criteria, considering both p-values and false discovery rate (FDR) q-values below the threshold of 0.05 to denote statistical significance. This approach ensured a robust and reliable interpretation of our data, minimizing the likelihood of false positives and providing a solid foundation for our conclusions.

The graphic flowchart in Fig. 1 succinctly encapsulates the methodology and structure of our extensive study on the role of SMRGs in HNSCC. Our research entailed a detailed analysis of clinical information and mRNA expression data from 546 HNSCC samples, which were sourced from The Cancer Genome Atlas (TCGA). In the initial phase of our study, we identified a comprehensive gene set of soluble mediators, totaling 10,859 genes, from the Genecard database. This extensive list served as a foundational reference for our subsequent analyses. Employing the “limma” R package, we conducted differential gene expression analysis on HNSCC tumor tissues. We set a stringent threshold of an absolute log${}_2$ fold change ($|$log${}_2$FC$|$) greater than 1, which led to the identification of 1943 upregulated and 537 downregulated genes (Fig. 2A). To visually represent these findings, a heatmap was generated, showcasing the mRNA expression differences between tumor samples (n = 502) and normal samples (n = 44) (Fig. 2B). This visualization provided an immediate and clear comparison of gene expression patterns between the two sample types. Further, we intersected the differentially expressed genes identified in tumor samples with the initial soluble mediator gene set. This intersection resulted in the identification of 2361 differentially expressed SMRGs (DE-SMRGs). These DE-SMRGs were then subjected to further analysis using the Clusterprofiler R package for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment. The criteria for significant enrichment were set at an adjusted p-value (p.adj) of less than 0.05 and a false discovery rate (FDR) of less than 0.25. The results of this enrichment analysis were presented in bubble maps (Fig. 2C,D), highlighting the top 25 GO and top 20 KEGG enrichment results. These results predominantly showed enrichment in pathways related to T cell activation, lymphocyte regulation, and immune responses, including Th1 and Th2 cell differentiation, T cell receptor signaling, TNF signaling, and IL-17 signaling pathways. The enrichment in these specific pathways underscores the potential immunological functions and implications of SMRGs in the pathology of HNSCC.

Fig. 1.

The workflow of identification of the SMRGs signature for patients with HNSCC. PCA, principal components analysis; TCGA, The Cancer Genome Atlas; SMRGs, soluble mediator-related genes; ICGC, International Cancer Genome Consortium; TISCH, Tumor Immune Single-Cell Hub; ROC, Receiver Operating Characteristic; HNSCC, Head and Neck Squamous Cell Carcinoma; K-M, Kaplan-Meier.

Fig. 2.

Biological process analysis of soluble mediators related genes. (A) Volcano map of DE-SMRGs. (B) DE-SMRGs in HNSCC tumor samples and normal samples. (C) GO enrichment analysis DE-SMRGs. (D) KEGG pathway analysis of DE-SMRGs. DE-SMRGs, differentially expressed SMRGs; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change.

We categorized the TCGA-HNSC cohort into groups based on expression of DE-SMRGs using the ConsensusClusterPlus R program (Supplementary Figs. 1,2). A consensus matrix with k value equal to 2 was found to minimize crossover between HNSC samples (Fig. 3A). Survival differences between clusters were analyzed using the ClusterSurvival R package, revealing better overall survival (OS) in one cluster compared to the other (Fig. 3B). To elucidate the reasons behind these prognostic differences, we examined the levels of immune-related cells in both clusters using the MCPcounter R package, including B lineage, CD8 T cells, cytotoxic lymphocytes, endothelial cells, monocyte lineage, bone marrow dendritic cells, NK cells, and T cells. We found significantly higher infiltration of these cells in one cluster (Fig. 3C–J), suggesting a correlation between immune cell infiltration and better patient prognosis in HNSCC.

Fig. 3.

Immunotyping of differentially expressed genes associated with soluble mediators. (A) Consensus clustering matrix for K = 2. (B) Survival curves of different immune genotypes. (C–J) The content of important immune cells among different types.

In our study, we developed a sophisticated risk score model based on SMRGs to predict outcomes in HNSCC patients. Initially, univariate Cox analysis, executed via the “survival” R package, identified 274 prognostic-related SMRGs (p 0.05). Further refinement through lasso regression analysis distilled these to 45 key SMRGs (Fig. 4A,B, Supplementary Table 1). Employing a multifactorial Cox proportional hazard regression model, we narrowed down to 20 pivotal SMRGs, specifically CSF2, ITGA4, LGALS3BP, EPO, EZH2, SFRP1, PSMA1, TNFRSF4, GRIA3, CYP2D6, SLC25A4, CCNA1, PRPS1, UCN2, CCL28, HEXA, SPINK1, HTN3, TLL1, and SPINK6 (Fig. 4C). The corresponding regression coefficients for these genes were calculated as follows: 0.3415, –0.5791, 0.9097, –0.7959, –0.5625, –0.1741, 0.8485, –0.8744, –0.8864, –0.7582, 0.5207, 0.3389, 1.4621, –0.3439, –0.3837, 1.2022, 0.4666, 0.484, 0.4744, and –0.3499. This linear prediction model, derived from multivariate Cox analysis, based on 20 weighted SMRGs, allowed us to calculate a patient-specific risk score.

Fig. 4.

Construction and validation of SMRGs Signature and its prognostic value. (A) Ten-time cross-validation for tuning parameter selection in the LASSO model. (B) LASSO coefficient profiles. (C) Prognosis of 20 SMGs in the entire cohort of HNSC was analyzed by univariate Cox regression model. K-M survival curve of HNSC patients in the TCGA train cohort (D), TCGA test cohort (E) and ICGC test cohort (F). Time-dependent ROC curve of HSNC patients in the TCGA train cohort (G), TCGA test cohort (H) and ICGC test cohort (I). (J) Genetic correlation analysis of 20 SMGs. (K,L) Pathway analysis of 20 SMGs. SMGs, soluble mediator-related genes; HNSCC, Head and Neck Squamous Cell Carcinoma.

To ensure the robustness of our prognostic model, we divided the patients into a training cohort (n = 353) and a testing cohort (n = 148), following a 7:3 ratio. The training cohort results aligned with the overall findings, where HR patients had significantly worse outcomes (p 0.001) (Fig. 4D), and the ROC curve area under the curve (AUC) values for 1, 3, and 5-year predictions were 0.745, 0.838, and 0.839, respectively (Fig. 4G). This trend was consistent in the testing cohort, with HR patients demonstrating poorer prognoses (p = 0.017) (Fig. 4E), although the ROC curve AUCs at 1, 3, and 5 years were 0.708, 0.635, and 0.564, respectively (Fig. 4H). For external validation, we utilized 40 samples from the International Cancer Genome Consortium (ICGC) database as an independent test cohort. The results were consistent with the earlier findings, indicating a poorer prognosis for HR patients (p = 0.022) (Fig. 4F), and the ROC curves showed AUCs of 0.850 at 1 year and 0.870 at 1.5 years (Fig. 4I).

To explore the biological mechanisms behind the aberrant expression of the 20 target genes, we conducted gene correlation analysis (Fig. 4J) and pathway analysis. The pathway analysis revealed significant involvement of these SMRGs in regulating critical processes such as apoptosis, epithelial-mesenchymal transition (EMT), DNA damage response, and the cell cycle, among other tumor-related pathways (Fig. 4K). We further classified these pathways into activated (A) and repressed (I) categories, as depicted in Fig. 4L. Notably, these genes were found to influence the activation of the PI3K/AKT pathway in HNSCC patients. Other pathways, including apoptosis, EMT, RTK, RAS/MAPK, cell cycle, Hormone ER, Hormone AR, and DNA damage response, exhibited varying degrees of activation or repression.

We undertook a detailed clinical subgroup analysis to evaluate the variations in patient prognosis across diverse clinical parameters. This analysis was essential to understand how different patient characteristics might influence outcomes when assessed through the lens of our SMRGs risk model. To this end, we meticulously stratified the HNSCC patient cohort into various subgroups. These subgroups were based on critical clinical criteria such as T stage (T1-2 and T3-4), pathological stage (I–III and III–IV), pathological N stage (grades N0 and N1–3), tumor grade (grades I–II and III–IV), age ($\ge$65 years), and gender (male and female), as illustrated in Fig. 5. The stratification allowed us to conduct a nuanced and comprehensive analysis of survival outcomes across these different subgroups. Remarkably, the results consistently indicated that patients in the high-risk (HR) group, as determined by our SMRGs risk model, had significantly poorer overall survival (OS) in all the clinical subgroups examined (Fig. 5). This finding not only validates the predictive strength of our risk model but also highlights its applicability and robustness across various clinical scenarios in HNSCC. It underscores the potential of the SMRGs-based risk model as a powerful prognostic tool capable of guiding treatment decisions and predicting outcomes in a wide range of clinical contexts within HNSCC management.

Fig. 5.

Prognostic power of the SMRGs risk model for overall survival for multiple HNSC subtypes. (A) Age $>$65 years. (B) Age $\le$65 years. (C) Female. (D) Male. (E) Grade 1–2. (F) Grade 3–4. (G) N0. (H) N1–3. (I) Stage I–II. (J) Stage III–IV. (K) T1–2. (L) T3–4.

Our in-depth study on HNSCC entailed rigorous univariate and multivariate Cox regression analyses to evaluate the prognostic potential of our SMRGs-based model. These analyses were crucial for determining whether our 20-SMRGs-based prognosis signature could act as an independent prognostic factor, considering the overall survival (OS) of HNSCC patients in conjunction with their clinical characteristics. The univariate analysis yielded significant findings, revealing a strong correlation between several clinical parameters—age, stage, grade, T, N—and the risk score with patient prognosis (p 0.001) (Fig. 6A). This association underscored the relevance of these clinical factors in determining patient outcomes. Moving to the multivariate analysis, the risk score continued to stand out as a robust and independent prognostic indicator within our patient cohort (p 0.01) (Fig. 6B). To further illustrate these relationships, a heatmap was generated, showcasing the association between the expression levels of the 20 identified SMRGs in our prognostic risk model and various clinical parameters—age, grade, stage, T, N, and risk score—across all HNSCC patient samples from the TCGA (Fig. 6C). In an effort to enhance the clinical utility of our risk model, we developed a comprehensive nomogram. This nomogram incorporated a range of variables: age, grade, stage, T, N, the expression levels of the 20 SMRGs, and the risk score, to predict 1, 3, and 5-year survival probabilities for HNSCC patients (Fig. 6D). The results from the nomogram clearly highlighted the dominance of the risk score in predicting OS, suggesting that our 20-SMRGs-based model offers a more precise prognostic prediction for HNSCC compared to conventional clinical parameters. The accuracy and reliability of this nomogram were further validated through calibration curves, which showed a close match between the predicted and observed values for 1-year, 3-year, and 5-year OS (Fig. 6E). Notably, the AUC value of 0.731 for the risk score derived from the 20 SMRGs demonstrated its superior predictive power for HNSCC prognosis, surpassing that of common clinicopathological features such as age, grade, pathological stage, T, and N (Fig. 6F). This finding reinforces the model’s significance in clinical decision-making, evidenced by its high net benefit depicted in the decision curve analysis (Fig. 6G).

Fig. 6.

Independent prognostic analysis of risk scores and clinical parameters. Univariate (A) and multivariate (B) COX regression analysis of the signature and different clinical features. (C) Heatmap for the 20 SMRGs-based signatures with clinicopathological manifestations. (D) Nomogram for predicting 1-year, 3-year, and 5-year OS of patients with HNSC. (E) The calibration curve of the constructed nomogram of 1-year, 3-year, and 5-year survival. (F) Multi-index ROC analysis in the test cohort. (G) Decision curve analysis. OS, overall survival; AUC, area under the curve.

We developed a novel SMRGs signature and compared its predictive power against four established HNSCC prognostic signatures: the Lian signature [51], Tang signature [52], Xue signature [53], and Huang signature [54] (Fig. 7A–D). Uniform risk scores were calculated for each HNSCC sample across all TCGA cohorts to ensure a standardized comparison. Time-dependent ROC curve analysis revealed that our SMRGs signature consistently outperformed the other signatures, exhibiting higher AUC values at 1-year, 3-year, and 5-year survival intervals. Notably, it achieved the highest C-index of 0.718 (Fig. 7E), demonstrating superior predictive performance. This underscores the potential of our SMRGs signature as a more effective prognostic tool in HNSCC, offering enhanced accuracy in predicting patient outcomes and informing personalized treatment strategies.

Fig. 7.

Comparison of the SMRGs signature with other models. (A–D) Kaplan-Meier (K-M) curves and ROCs for risk models constructed by others. (E) C-indexes for five risk models.

We explored the correlation between our SMRGs-based risk models and key clinicopathological characteristics of patients. Using the Wilcoxon test, we analyzed variations in risk scores across subgroups defined by gender, T stage, pathological stage, tumor grade, and age. The results indicated no significant association between risk scores and age, N stage, or gender. However, a notable correlation was found with tumor grade (p 0.01), T stage (p 0.05), and pathological stage (p 0.05) (Fig. 8). This analysis highlights the relevance of our SMRGs-based risk model, particularly in relation to tumor severity and progression in HNSCC, offering valuable insights for patient management and treatment planning.

Fig. 8.

Correlation analysis of risk scores and clinicopathological characteristics. (A) Age. (B) Gender. (C) Tumor grade. (D) Stage. (E) T stage. (F) N stage.

We conducted an in-depth investigation into the relationship between the SMRGs risk score and the tumor microenvironment (TME) in HNSCC, focusing particularly on immune cell infiltration. Utilizing a comprehensive set of algorithms, including XCELL, QUANTISEQ, TIMER, MCPCOUNTER, CIBERSORT-ABS, CIBERSORT, and EPIC (Fig. 9A), we sought to understand the complex dynamics within the TME. Our findings indicated a significant inverse relationship between the SMRGs risk score and the presence of various immune cells, including myeloid dendritic cells, T cells, CD8 T cells, B lineage cells, and NK cells (Fig. 9B). Moreover, single-sample Gene Set Enrichment Analysis (ssGSEA) revealed pronounced differences in immune function between high-risk (HR) and low-risk (LR) groups (Fig. 9C). A particularly striking discovery was the differential expression of immune checkpoints between the HR and LR groups. This aspect is of critical importance in the context of checkpoint-based immunotherapy. The LR group showed heightened expression of twelve immune checkpoint genes, notably including LAIR1, IDO1, CTLA-4, PD-1, TIGIT, CD200R1, CEACAM1, CD200, KIR3DL, BTLA, and ADORA2A, with a pronounced expression of PD-1. This pattern suggests a dependency on the PD-1/PD-L1 pathway for immune evasion in LR tumor cells (Fig. 9D), implying that treatments targeting PD-1 might be particularly effective in this subgroup. The activation of these immune checkpoints, such as PD-1, is typically associated with an inflamed TME [55], suggesting an inflammatory milieu in LR patients. Complementing this, our comparison of ssGSEA scores for various immune cells and functions revealed that all immune cell types, except for macrophages, were significantly more abundant in the LR group. Additionally, nine immune functions, including HLA and CCR, scored higher in the LR group (Fig. 9E,F). This comprehensive analysis underscores not only the robustness of our SMRGs risk score in prognostic prediction for HNSCC but also its significant implications in shaping personalized immunotherapeutic strategies. The insights garnered from this study could potentially guide more effective and tailored treatments for patients with HNSCC.

Fig. 9.

SMRGs risk score predicts tumor microenvironment and immune cell infiltration. (A) Immune cell bubble plots. (B) The risk score and immune cell correlation. (C) Immune cell ssGSEA scores between high and low-risk groups. (D) Immune checkpoint differences between high and low-risk groups. (E) Immune function ssGSEA scores between high and low-risk groups. (F) Heat map of immune function differences between high- and low-risk groups. *p 0.05; **p 0.01; ***p 0.001. ssGSEA, single-sample Gene Set Enrichment Analysis; ns, not significant.

We initiated our analysis by exploring the Immunotherapy Prediction Pathway, focusing on the comparison between high-risk (HR) and low-risk (LR) groups. We discovered significant differences in several pathways (Fig. 10A). Notably, the LR group exhibited elevated Systemic lupus erythematosus risk scores, while other pathways such as Homologous recombination, Spliceosome, Progesterone-mediated oocyte maturation, and RNA degradation presented lower risk scores compared to the HR group (Fig. 10B). This suggests a better prognosis in the LR group, potentially influenced by aberrant RNA degradation processes impacting gene expression and cancer progression [56]. Further, we utilized the tumor immune cycle as a framework to evaluate the impact of immunomodulators and chemokine systems [44, 57]. Our comparison of the HR and LR groups showed marked differences in the activities of various steps of the tumor immune cycle. Activities related to the expression of cancer cell antigen (step 2), initiation and activation (step 3), transport of immune cells to the tumor (step 4), infiltration into tumors (Step 5), recognition of cancer cells by T cells (Step 6), and the killing of cancer cells (step 7) were predominantly upregulated in the LR group. In contrast, the release of cancer cell antigens (step 1) demonstrated decreased activity (Fig. 10C). Analysis of risk scores related to these differences indicated higher scores in steps 3, 4, 5, 6, and 7 in the LR group, while other steps showed no significant variation (Fig. 10D).

Fig. 10.

SMRGs risk scores predicting treatment response assessment. (A) Heat map of the difference in enrichment scores between high and low-risk groups on the immunotherapy prediction pathway. (B) Heat map of differences between high and low-risk groups on each step of the cancer-immune cycle. (C) Differences in risk scores of immunotherapy prediction pathway between high and low-risk groups. (D) Differences in risk scores between high- and low-risk groups at steps of the cancer-immune cycle. (E) Correlation of different risk scores with the corresponding signals at the cancer-immune cycle step and ICB. (F–J) Sensitivity analysis of SMRGs score in predicting immunotherapy response to cancer. *p 0.05; **p 0.01; ***p 0.001. ICB, immune checkpoint blockade.

Additionally, we explored the correlation between immune checkpoint blockade (ICB)-related signals and SMRG risk scores. This analysis revealed a negative correlation of the risk scores with pathways including DNA replication, Cell cycle, Fanconi anemia pathway, Homologous recombination, Nucleotide excision repair, Mismatch repair, Oocyte meiosis, Progesterone-mediated oocyte maturation, Spliceosome, Pyrimidine metabolism, and Systemic lupus erythematosus, while a positive correlation was observed with RNA degradation. Other pathways did not show significant correlations (Fig. 10E). Moreover, we found that SMRG risk scores were inversely related to most tumor immune processes, except for cancer cell antigen release, Th22, eosinophil, neutrophil, basophil, or MDSC cell recruitment (Fig. 10E). Finally, our analysis using the BEST database demonstrated that the SMRGs Score has high sensitivity in predicting responses to immunotherapy in cancer patients (Fig. 10F–J). This finding provides valuable insights into the potential of SMRGs as a tool for guiding immunotherapeutic strategies, underscoring the significance of our comprehensive approach in understanding and predicting HNSCC treatment outcomes.

A.770041 (p = 3.4 $\times$ 10${}^{-13}$), AS601245 (p = 6.9 $\times$ 10${}^{-10}$), AZ628 (p = 4.8 $\times$ 10${}^{-11}$), and AZD.0530 (p = 1.9 $\times$ 10${}^{-5}$) had greater IC50s in the LR group compared to the HR group when used to treat HNSC, out of a total of nine immunotherapeutic agents (Fig. 11A,F,H,I). The IC50 for the other 5 chemical or targeted medicines was lower in the LR group: ABT.263 (p = 0.024), ABT.888 (p = 2.1 $\times$ 10${}^{-7}$), AICAR (p = 5.4 $\times$ 10${}^{-6}$), AKT inhibitor VIII (p = 1.9 $\times$ 10${}^{-7}$), and ATRA (p = 9.1 $\times$ 10${}^{-11}$) (Fig. 11B,C,D,E,G). The risk score allows for more in-depth research into the immunotherapy response in HNSC patients and more targeted medication therapy.

Fig. 11.

Differences in IC50 of immunotherapy drugs by risk score. (A) A.770041. (B) ABT.263. (C) ABT.888. (D) AICAR. (E) AKT. inhibitor Ⅷ. (F) AS601245. (G) ATRA. (H) AZ628. (I) AZD.0530. IC50, half-maximal inhibitory concentration.

We leveraged the single-cell dataset HNSC GSE103322 from the TISCH database to explore the expression profiles of 20 selected SMRGs in tumor microenvironment (TME)-related cells. The dataset encompasses 20 distinct cell clusters, representing 11 unique cell types, as detailed in Fig. 12A. Our analysis revealed a varied expression pattern of SMRGs among different cell types. Notably, a spectrum of immune cells, including proliferating T, myofibroblasts, CD8${}^{+}$ T, CD4${}^{+}$ T, and B cells, exhibited pronounced expression of LGALS3BP, PSMA1, HEXA, and PRPS1. Conversely, malignant cells predominantly expressed SPINK6, TLL1, UCN2, CCNA1, EPO, and CCL28. Proliferating T cells were characterized by elevated levels of ITGA4, TNFRSF4, PRPS1, and CYP2D6. Additionally, we observed a heterogeneous expression of PSMA1 and EZH2 across various immune cell types. LGALS3BP and HEXA, while displaying low to moderate expression in several immune cell populations, were markedly expressed in fibroblasts, myofibroblasts, and malignant cells. Furthermore, SFRP1 was notably present in fibroblast cells, whereas SLC25A4 was primarily found in myofibroblasts. Lastly, fibroblasts and cancer cells exhibited substantial expression of CSF2, GRIA3, HTN3, and SPINK1 (Fig. 12B, Supplementary Fig. 3).

Fig. 12.

SMRGs Expression in HNSCC TME-associated cells. (A) Annotation of all cell types in GSE103322 and the percentage of each cell type. (B) Percentages and expression of 20 SMRGs in different cell types in GSE103322. TME, tumor microenvironment.

To further validate gene expression, we used immunohistochemical analysis. The results showed that LGALS3BP, EZH2, SFRP1, CCNA1 and SPRINK1 were highly expressed in HNSCC, while SLC25A4 was highly expressed in normal samples (Fig. 13).

Fig. 13.

Immunohistochemical analysis was used to verify gene expression. (A) LGALS3BP. (B) EZH2. (C) SFRP1. (D) SLC25A4. (E) CCNA1. (F) SPRINK1. Images available from https://www.proteinatlas.org/ENSG00000108679-LGALS3BP/pathology/head+and+neck+cancer#img. Light yellow, weak positivity; brownish yellow, moderate positivity; brownish black, strong positivity