GSEA was performed using the GSEA software (version 4.2.3) (https://www.gsea-msigdb.org/gsea/downloads.jsp) with the MSigDB glycolysis-related pathway gene sets, which contain 1320 gene sets. Pathways with p < 0.05 and FDR < 0.05 were considered significantly enriched.

The HNSCC gene expression data (RNA-Seq) and the corresponding clinical data (including age, gender, stage, grade, smoking, alcohol, HPV, survival time, and survival status) were downloaded from the TCGA database (https://portal.gdc.cancer.gov) and GEO dataset (https://www.ncbi.nlm.nih.gov/geo/). Used as a training cohort, the inclusion criteria for TCGA-HNSCC were: HNSCC samples with complete somatic mutation data and clinical information (457 retrieved HNSCC samples with single nucleotide polymorphism(SNP) data were analyzed), with 462 HNSCC samples and 32 adjacent non-tumor tissue samples included. The glycolysis-related gene dataset was downloaded from MSigDB. Expression characteristics of glycolysis-related genes were obtained from the MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).

Using |log FC| > 1 and p < 0.05 as thresholds, differentially expressed genes between HNSCC samples and adjacent non-tumor tissue samples were evaluated using the Wilcoxon test in the limma package. Then, differentially expressed glycolysis-related genes were selected from all differentially expressed genes (DEGs) and displayed on a Venn diagram.

Survival-associated differentially expressed glycolytic genes were identified via univariate Cox regression and Lasso regression analyses, after which a polygenic prognostic risk model was constructed. Based on the median risk score of the TCGA training set as the cutoff, HNSCC patients were assigned into high- and low-risk groups. Clustering effects of Principal Component Analysis (PCA) dimensionality reduction revealed significant differences between the groups. Kaplan-Meier survival curves, time-dependent receiver operating characteristic (ROC) curves, and risk score distributions for OS prediction were evaluated to verify the prognostic significance of risk scores. Similar to the training set approach, the GEO cohort was used as an independent validation set to assess the generality and reliability of the prognostic risk model.

Independent prognostic factors in HNSCC patients were determined by univariate and multivariate Cox regression analysis. Both TCGA training set and GEO validation set were used to construct a nomogram for predicting the 1-year, 3-year, and 5-year survival outcomes of HNSCC patients. Consistency between actual survival rates and nomogram-predicted rates was tested via a calibration curve. In addition, decision curves were used to assess the reliability of risk scores and clinical stage.

The “Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT)”was used to assess immune cell infiltrations. The immune, stromal, and ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumors using Expression data) scores for each sample were calculated using the ESTIMATE algorithm. Correlations between the GRGPI score and those scores were determined by Spearman correlation analysis.

The tumor mutation data was obtained from the cBioPortal database. The tumor mutation burden (TMB) for all samples was calculated using “maftools” in R. Based on median TMB values, HNSCC samples were assigned into high TMB and low TMB groups. A total of 16360 genes involved in developing SNP in 457 samples were obtained by MusigCV (running under the linux system), and the top ten were screened using q<0.05 as the cut-off. Correlations between the prognosis for HNSCC patients with GRGPI and TMB were determined by Kaplan-Meier survival curves in R.

To assess the clinical applicability of the established model, pRRophetic was used to calculate the IC50 of HNSCC chemotherapeutic drugs.

The R software (version 4.1.1) was used for statistical analyses. Differentially expressed genes between tumor and adjacent normal tissues were compared by the Wilcoxon test. Survival-associated differentially expressed glycolytic genes were identified by univariate Cox and Lasso regression analyses. Then, Kaplan-Meier survival curves were plotted. Univariate COX and multivariate COX regression analyses were performed to determine the independent prognostic factors for OS. The predictive ability of the model was assessed by KM survival curves and ROC curves. Correlation tests were conducted by Spearman correlation analyses. Categorical data were compared by the chi-square test. p ≤ 0.05 was the threshold for statistical significance. The flow chart of our study is shown in Figure 1 .

All HNSCC and adjacent non-tumorous tissue samples were collected from 10 patients in the Shanxi Province Cancer Hospital. Extraction of total RNA from HNSCC tissues and adjacent non-tumor tissues was performed by the TRIzol reagent (Invitrogen, CA, USA). cDNA synthesis from the extracted RNA was performed by PrimeScriptTM RT Master Mix (RR036B, Takara). We use Quantitative PCR to analyze the mRNA expression levels of the prognostic glycolysis-related genes by GoTaq^® qPCR Master Mix (Promega, A6001). The RT-qPCR was utilized in the ABI Vii7 Sequence detection system (ABI, USA). The PCR reaction system and conditions were according to the manufacturer’s instructions. Gene expression levels of STC1, STC2, AURKA, P4HA1, and PLOD2 were calculated using the 2-ΔΔCT method.

Based on KEGG, REACTOME and HALLMARK gene sets, GSEA was performed to reveal potential differences between HNSCC and control groups. These pathways are associated with glycolysis, implying that glycolysis plays an essential role in HNSCC ( Supplementary Figures 1A–C ).

A total of 1695 differentially expressed genes (DEGs) (149 upregulated and 119 downregulated genes) in the TCGA training cohort were identified by Wilcoxon signed-rank test and visualized using volcano plots ( Figure 2A ) and heatmaps ( Figure 2B ). Then, 49 glycolysis-related genes were extracted from the DEGs ( Figure 2C ).

Through univariate Cox regression analysis, 15 prognosis glycolytic genes were established to be closely associated with survival outcomes of HNSCC patients ( Figure 3A ). The 15 OS-related genes may be collinear rather than independent. LASSO Cox regression analysis was performed to determine the real OS-affecting factors, and finally, a prognostic panel consisting of five glycolysis-related genes was established. The risk score was calculated as: Riskscore=0.021*AURKA+0.099*P4HA1+0.015*PLOD2+0.031*STC1+0.163*STC2. ( Figure 3B, C ). Based on this gene signature, all patients were assigned to high (n=231) and low-risk (n=231) subgroups using the risk score median as the threshold. Risk scores, survival scores, and heatmap of prognostic glycolytic gene expressions among the low-risk and high-risk patients are presented in Figure 3D . Based on expressions of these five hub genes, dimensionality reduction was performed in all patients and presented with methods of t-distributed stochastic neighbor embedding (t-SNE), suggesting that different risk subgroups show significant discrete tendencies directly in the two-dimensional plane ( Figure 3E ). Kaplan-Meier survival curves revealed that high-risk score patients had significantly worse OS outcomes than low-risk score patients. The area under the curve (AUC) analysis for HNSCC patients at 1-year, 3-year, and 5-year revealed respective OS rates of 0.622, 0.649, and 0.614, demonstrating the optimal predictive performance of GRGPI ( Figures 3F, G ). Finally, the year with the largest AUC value is shown in the RMST plot ( Figure 3H ).

Given that the predictive potential of GRGPI in different datasets is misty into account, GSE65858 was used as the independent validation set. Based on the above risk scores, patients were assigned to low-risk (n=140) and high-risk (n=130) groups ( Figure 4A ). Findings from t-SNE and KM survival analyses of the GEO validation set were consistent with the results of the TCGA training cohort ( Figures 4B, C ). The AUC values for ROC curves accurately revealed the predictive performance of the prognostic risk model, with the largest AUC value year shown as an RMST plot ( Figures 4D, E ).

In this study, the risk score, gender, smoking, and clinical stage were established to be independent prognostic factors in HNSCC patients, and they were used to construct subsequent nomograms ( Figure 5A ). Nomograms were used to predict the 1-year, 3-year, and 5-year survival probabilities of HNSCC patients ( Figure 5B ). Moreover, a calibration curve was constructed to assess the agreement between nomogram predictions and actual survival outcomes ( Figure 5C ). The actual and predicted survival rates at 1-year, 3-year, and 5-year were well matched, indicating that the nomogram has a good predictive performance. A decision curve (DCA) was used to assess the reliability of the risk score. It was observed that Model1 (Stage) was close to the extreme curve, while Model2 (RiskScore) was significantly higher than the extreme curve ( Figure 5D ).

The above analyses were also performed on the GEO validation set to verify the robustness of the model ( Figure 6A ). Unlike the TCGA training set, univariate and multivariate Cox regression analyses revealed that in the GEO validation set, only risk score and clinical stage were independent prognostic factors for HNSCC patients. Therefore, a nomogram integrating risk scores and clinical stages was constructed to predict the 1-year, 3-year, and 5-year survival probabilities of HNSCC patients ( Figure 6B ). Findings from the calibration curve were consistent with those of the TCGA training set ( Figure 6C ).

Based on the relationship between high and low-risk groups and clinical stages in the TCGA training set, a clinical correlation table was prepared. It was established that about 60% of patients in the low-risk group were in locations I/II, while 76% of patients in the high-risk group were in stages III/IV ( Supplementary Figure 2 ), implying that risk grouping was positively correlated with the clinical stage. These findings prove that the constructed GRGPI model is clinically valuable.

Mutation data for HNSCC were downloaded from the cBioPortal database. Somatic mutation types for the 457 patients were evaluated, and SNP was found to be the most dominant mutation type ( Figures 7A, C ). About 95.18% of samples in the high-risk group had SNPs, compared to 96.07% in the low-risk group ( Figures 7B, D ).

Given that SNP was the most dominant mutation type in HNSCC, 16,360 SNP-mediating genes were identified in the 457 samples using MusigCV. Using q<0.05 as the cut-off, the top ten genes were screened. Mutations types of the 10 genes and their distributions in high-risk and low-risk groups were analyzed ( Figure 7E ). The top ten genes with the highest mutation rates in the high-risk group were TP53, TTN, FAT1, MUC16, CSMD3, CDKN2A, LRP1B, KMT2D, DNAH5, and PIK3CA ( Figure 7B ), while those in the low-risk group were TP53, TTN, FAT1, MUC16, CSMD3, NOTCH1, SYNE1, CDKN2A, LRP1B, and PIK3CA ( Figure 7D ). It was observed that the gene with the highest mutation rate was TP53 in HNSCC patients regardless of the high GRGPI group or the low GRGPI group, suggesting that the mutations of the tumor suppressor gene TP53 may have potential clinical and pathophysiological significance in HNSCC patients. In fact, in a recent study, the mutational profile of TP53 has been proved to act as an independent prognostic factor in HNSCC patients. This relationship is associated with unique site-specific biological networks, consistent with our findings (17). Correlation analyses showed that GRGPI was positively correlated with TMB (R=0.015, p=0.75). The difference in the number of HNSCC patients in the high and low TMB groups was insignificant ( Supplementary Figure 3A ). Moreover, the difference in TMB values between the groups was negligible ( Supplementary Figure 3B ). We combined GRGPI and TMB and grouped them into three; high GRGPI high mutation (HTMB+HGRGPI), high GRGPI low mutation or low GRGPI high mutation (HTMB+LGRGPI & LTMB+HGRGPI), and low GRGPI low mutation (LTMB+LGRGPI). Then, Kaplan-Meier survival curves were drawn. The survival curve showed that the LTMB+LGRGPI group had the best prognosis, while the HTMB+HGRGPI group had the worst prognosis ( Figure 7F ). These findings imply that high GRGPI and high TMB play a synergistic role in promoting tumor occurrence and development, and the combined effects of the two may lead to worse prognostic outcomes.

Interactions among the five glycolysis key genes and transcription factors may elucidate on mechanisms of the GRGPI model. Using cor Filter=0.5 and fdr Filter=0.01 as critical values, associations between AURKA, P4HA1, PL0D2 and 11 transcription factors were obtained ( Supplementary Figure 4 ), which proved that the genes used to construct the GRGPI model were correlated with transcription factors in cancer and para cancer differentially expressed genes.

Compared to the low-risk group, infiltrations of resting CD4 memory T cells, M0 macrophages, M2 macrophages, and activated mast cells were marked in the high-risk group, while infiltrations of CD8 T cells, follicular helper T cells, and Treg cells were to a greater extent ( Figure 8A ). Cell immunity-related cells, such as CD8 T cells, were highly infiltrated in the low-risk group, suggesting that immune cells may be activated in the low-risk group and suppressed in the high-risk group. Moreover, the M0 macrophages, activated mast cells, and resting CD4 memory cells were positively correlated with GRGPI scores while resting dendritic cells, CD8 T cells, follicular helper T cells, and Treg cells were negatively correlated with GRGPI scores ( Figure 8B ). The higher the GRGPI scores, the worse the extent of T cell infiltrations, validating that weaker antitumor immunity may be one of the reasons for poor prognostic outcomes. Therefore, the high-risk group was defined as the low-immunity group, while the low-risk group was defined as the high-immunity group. Differences in ESTIMATE scores between the high-risk group and low-risk group were insignificant. However, the high-risk group exhibited low immune scores (p< 0.001, Figure 8D ). These findings are consistent with those obtained from CIBERSORT, whereby the high-risk group exhibited worse immune status while the low-risk group exhibited better immune status. The relationship between stromal cells and GRGPI scores was further investigated (18). The high-risk group had higher stromal scores (p < 0.01, Figure 8E ), implying that tumor stroma plays an important role in tumor development.

ICI therapy has advanced the treatment of many solid tumors. Therefore, 11 human leukocyte antigen(HLA) class immune checkpoints were included, and their differential expressions in high-risk and low-risk groups were determined. Four HLA class checkpoints (HLA-A, HLA-C, MICA, and MICB) were highly expressed in the high-risk group ( Figure 9A ), while the remaining seven were highly said in the low-risk group. Since the HLA class immune checkpoints are closely associated with immune responses, the better prognostic outcomes in the low-risk group could have been due to better immune responses. The expressions of 7 genes (CD274, CTLA4, IDO1 LAG3, PDCD1, TIGIT, and TNFRSF9) in the high-risk group and low-risk group were also analyzed ( Figure 9B ). Results show that five immune checkpoints cut in the high-risk group, CTLA4, IDO1, LAG3, PDCD1, and TIGIT, which are consistent with the result, once again proved that the GRGPI model and the close correlation between HLA class immune checkpoints.

GSEA was performed to assess the immune pathways, and differentially expressed immune-related pathways between the high-risk and low-risk groups were obtained ( Figure 9C ). The BCR, Chemokines, Chemokine, Receptors, Interleukins Receptor, NK cell Cytotoxicity and TCR signalling pathway were found to be enriched in the high-risk group. However, enrichments of “TGFβ Family Members” and “TGFβ Family Members Receptor” were significantly high in the low-risk group, in accordance with the functions of TGF-β, which is involved in tumorigenesis and immunosuppression.

The Tumor Inflammation Signature (TIS) is investigational use only (IUO) 18-gene signature that measures pre-existing but suppressed adaptive immune responses within tumors (19). The high-risk group had a low TIS score, implying that this group had weaker adaptive immune responses and worse prognostic outcomes ( Figure 9D ).

The CYT score is a novel cancer immune index calculated from mRNA expressions of GZMA and PRF1 (20). The transcriptional levels of GZMA and PRF1 were determined to assess the cytolytic activities of immune lymphocytes in HNSCC. Based on previous risk grouping, the low-risk group exhibited a higher CYT score ( Figure 9E ), implying that immune cells in the low-risk group had stronger cytolytic activities and anti-tumor immune response, leading to a better prognosis.

TIDE is a computational framework developed by Peng Jiang et al. to identify two tumor immune escape mechanisms (21). A higher TIDE score means a greater likelihood of immune evasion, indicating that a patient is less likely to benefit from ICI therapy and may have worse prognostic outcomes. The TIDE website was used to process 457 HNSCC samples with complete somatic mutation data in the training cohort, of which 131 responded to immunotherapy while the remaining 326 did not. Then, the GRGPIs of responding and non-responding samples were evaluated, which revealed that responding samples had lower GRGPIs ( Figure 10A ). This confirms our findings in a previous study. Since the low-risk group had better performance in immune gene expressions, immune cell infiltrations, and activation of immune pathways, the higher degree of immune cell infiltrations enables it to achieve better results in immunotherapy, proving that our definition of the low-risk group as the high-immunity group in terms of immunotherapeutic effects is correct. Then, TMB values of response and non-response samples were determined, which did not reveal significant differences in TMB values ( Figure 10B ). The GRGPI established in this paper is superior to TMB in predicting immunotherapeutic effects. To validate the effects of immunotherapy in IMvigor210, differences in GRGPIs between response and non-response samples were investigated and found to be insignificant ( Figure 10C ). Differences between the two groups of TMB values were analyzed, and the response group was found to have higher TMB values ( Figure 10D ).

Therefore, the proportion of response and non-response samples in the three subgroups was identified after combining GRGPI and TMB ( Figure 10E ). The HTMB+HGRGPI group had the most significant proportion of responding models, followed by HTMB+LGRGPI& LTMB+HGRGPI, and LTMB+LGRGPI, suggesting that immunotherapy had better effects in the HTMB+HGRGPI group. To determine the prognostic performance of the established three subgroup models, we compared the AUC values of the three predictive models of GRGPI, TMB, and GRGPI combined with TMB ( Figure 10F ), which were 0.534, 0.647, and 0.646. The TMB and GRGPI combined with TMB exhibited better predictive performance. Finally, the prognostic value of the predictive model in melanoma was assessed using the GSE78220 cohort for external validation. Differences in GRGPIs between response and non-response groups were insignificant ( Supplementary Figure 5 ).

Although ICI therapy has shown great promise for the treatment of HNSCC, given its high costs and limited therapeutic effects (326/457 showed no responses to ICI therapy in this study), chemotherapy is a clinically meaningful treatment. However, HNSCC is associated with significant resistance to chemotherapeutic drugs during clinical treatment. To assess the application effects in the clinical chemotherapy process of the established model, IC50 was used to express the sensitivity of the high-risk and low-risk groups to several common chemotherapeutic drugs. Cisplatin, paclitaxel, and docetaxel were recommended by the Chinese Society of Clinical Oncology (CSCO) Guidelines of 2021 as first-line therapeutic drugs for HNSCC. Therefore, IC50 values in a high-risk group and low-risk group of the three drugs were calculated ( Figures 10G–I ). Patients in the high-risk group were more sensitive to cisplatin (p=1.4e-05) and docetaxel (p=5.5e-12). In contrast, those in the low-risk group were more sensitive to paclitaxel (p=9.9e-01), implying that the established model indicates chemotherapeutic sensitivity.

To verify the accuracy of GRGPI in HNSCC patients, we collected HNSCC tissues and adjacent non-tumorous tissues from 10 HNSCC patients. RT-qPCR was implemented to analyze the expressions of five prognostic glycolysis-related genes in the GRGPI. We found that compared with adjacent non-tumorous tissues, the terms of five genes were significantly up-regulated in HNSCC tissues ( Figures 11A–E ).