Skip to main content
Download PDF

Transcriptome analysis reveals the prognostic and immune infiltration characteristics of glycolysis and hypoxia in head and neck squamous cell carcinoma

BMC Cancer volume 22, Article number: 352 (2022) Cite this article

Abstract

Background

This study aims to construct a new prognostic gene signature in survival prediction and risk stratification for patients with Head and neck squamous cell carcinoma (HNSCC).

Method

The transcriptome profiling data and hallmark gene sets in the Molecular Signatures Database was used to explore the cancer hallmarks most relevant to the prognosis of HNSCC patients. Differential gene expression analysis, weighted gene co-expression network analysis, univariate COX regression analysis, random forest algorithm and multiple combinatorial screening were used to construct the prognostic gene signature. The predictive ability of gene signature was verified in the TCGA HNSCC cohort as the training set and the GEO HNSCC cohorts (GSE41613 and GSE42743) as the validation sets, respectively. Moreover, the correlations between risk scores and immune infiltration patterns, as well as risk scores and genomic changes were explored.

Results

A total of 3391 differentially expressed genes in HNSCC were screened. Glycolysis and hypoxia were screened as the main risk factors for OS in HNSCC. Using univariate Cox analysis, 97 prognostic candidates were identified (P < 0.05). Top 10 important genes were then screened out by random forest. Using multiple combinatorial screening, a combination with less genes and more significant P value was used to construct the prognostic gene signature (RNF144A, STC1, P4HA1, FMNL3, ANO1, BASP1, MME, PLEKHG2 and DKK1). Kaplan–Meier analysis showed that patients with higher risk scores had worse overall survival (p < 0.001). The ROC curve showed that the risk score had a good predictive efficiency (AUC > 0.66). Subsequently, the predictive ability of the risk score was verified in the validation sets. Moreover, the two-factor survival analysis combining the cancer hallmarks and risk scores suggested that HNSCC patients with the high hypoxia or glycolysis & high risk-score showed the worst prognosis. Besides, a nomogram based on the nine-gene signature was established for clinical practice. Furthermore, the risk score was significantly related to tumor immune infiltration profiles and genome changes.

Conclusion

This nine-gene signature associated with glycolysis and hypoxia can not only be used for prognosis prediction and risk stratification, but also may be a potential therapeutic target for patients with HNSCC.

Peer Review reports

Introduction

Head and neck squamous cell carcinoma (HNSCC), which includes a group of heterogeneous tumors from the squamous epithelium of the oral cavity, oropharynx, larynx and hypopharynx, is the seventh most common cancer in the world [1]. There are approximately 645,000 new cases of HNSCC each year worldwide. After receiving aggressive therapy, the 5-year survival rate of patients with HNSCC is still less than 50%. Patients with recurrent or metastatic HNSCC usually have a poor prognosis, with a median overall survival (OS) of only about 6 months. As far as we know, smoking, drinking, and human papilloma virus infection are risk factors for HNSCC. Moreover, the site of HPV infection is related to OS of HNSCC patients. For example, HPV positivity in the oropharynx, hypopharynx, oral cavity, and larynx is associated with improved OS, while HPV positivity in the nasopharynx and sinus passages is not associated with OS [2]. In the treatment of HNSCC, TNM staging is routinely used as the basis for selecting an appropriate treatment plan including surgery, radiotherapy and/or chemotherapy. Surgery is used for HNSCC patients whose tumors have not spread. Radiation therapy can be used for HNSCC patients with advanced pathological stage. Recent studies have shown that simultaneous radiotherapy and chemotherapy can improve patient survival, but it is not tolerated by many patients. Thus, although there have been some improvements in the treatment of HNSCC in the past few decades, the OS rate has not been improved significantly, which is mainly related to advanced stage at the time of diagnosis and the high treatment failure rate of advanced stage. Therefore, screening new prognostic biomarkers is essential to improve the clinical efficacy and OS of patients with HNSCC [3].

The hallmarks of cancer can show the basic characteristics of tumor cells. Recent studies suggested that some cancer hallmarks such as glycolysis and hypoxia were significantly related to the prognosis of patients with HNSCC. Glycolysis as the preferred pathway of energy metabolism is a characteristic of cancer cells. Therefore, cancer cells often exhibit increased glycolysis [4]. Previous studies have suggested that in HPV-negative HNSCC, the expression of genes related to glycolysis increased [5]. And glycolysis can be used as a biomarker to predict the prognosis of HNSCC patients [6]. In addition, pyruvate metabolism, a major glycolytic metabolite, was considered a potential anti-cancer target [7]. In terms of intratumoral hypoxia, tumor hypoxia is related to chemotherapy resistance and radiotherapy response in HNSCC [8, 9]. Moreover, hypoxia-inducible factor-1α was involved in intratumoral hypoxia-mediated tumor metastasis in HNSCC [10]. As hypoxic cells were more resistant to treatment, intratumoral hypoxia was an indicator of poor prognosis in HNSCC [11, 12]. Therefore, certain genes related to cancer hallmarks such as glycolysis and hypoxia are expected to serve as prognostic biomarkers for HNSCC.

Currently, treatment decisions for individual HNSCC patients are mainly based on the patient’s specific condition. However, the predictive ability and accuracy of traditional pathological staging have shown some deficiencies. Therefore, there is an urgent need to find new accurate and reliable predictors for HNSCC to guide risk stratification management and develop personalized treatment plans. In this study, we extracted an HNSCC cohort from the Cancer Genome Atlas (TCGA). Then, we explored the correlations between the cancer hallmarks and HNSCC. Next, we screened important prognostic genes related to glycolysis and hypoxia. A gene signature for survival prediction was then established. Subsequently, the prognostic value of the risk score based on gene signature was verified in the training set from TCGA and the validation set from the Gene Expression Omnibus (GEO). In addition, the correlations between risk scores and immune cell infiltration patterns, immune-related molecules and genomic changes were explored, respectively. Time-dependent receiver operating characteristic (t-ROC) curve was used to verify the prediction accuracy of the survival model. Taken together, this study comprehensively analyzed the prognostic value of a new gene signature related to cancer hallmarks including glycolysis and hypoxia in HNSCC. This gene signature can be used not only as a prognostic biomarker, but also as a potential therapeutic target for HNSCC.

Material and methods

Data set preparation and data processing

The training dataset and validation dataset for constructing prognostic gene signature were download from the TCGA and GEO databases, respectively. The training dataset with HNSCC-mRNA expression profile and clinical information obtained from the TCGA database (http://cancergenome.nih.gov/) included 546 HNSCC patients. The validation data sets (GSE41613 and GSE42743) with HNSCC-mRNA expression profile and clinical information downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) included 97 and 74 HNSCC patients, respectively. The above two databases are publicly available. Therefore, this study did not require the approval of the local ethics committee.

In order to make the gene expression profiles of different platforms comparable, we downloaded the data in FPKM format from the TCGA database and converted it into TPM format. Figure S1 showed the sample-normalized boxplots of two HNSCC cohorts from the GEO database. Subsequently, we performed a log2 normalization conversion on the above data. Meanwhile, using the limma package in R, the chip data downloaded from GEO was normalized [13]. Subsequently, 498 HNSCC samples with follow-up information and overall survival greater than zero from the TCGA RNAseq data were screened. The gene expression profile was obtained by removing the genes with zero expression levels in 50% of the samples. Meanwhile, the expression profile of the immune-related gene set was extracted. On the other hand, we performed the background correction on the gene chip data from GEO database. Then, 97 and 74 HNSCC samples with follow-up information and overall survival greater than zero from GSE41613 and GSE42743 were screened, respectively. Next, using the R package GEOquery, the chip probes were mapped to GeneSymbol. Finally, we get the gene expression profile by removing the probes mapped to multiple genes and taking the median of the multiple probes mapped to a single gene [13].

Selection of candidate genes and establishment of a gene signature

Based on the transcriptome profiling data and hallmark gene sets from the Molecular Signatures Database (MSigDB), the single-sample gene set enrichment analysis (ssGSEA) algorithm (R package “gsva”) was used to quantify the performance of cancer hallmarks in the training set [14, 15]. The R package “survival” was used to perform univariate Cox proportional hazards regression analysis (Cox-PH) to assess the significance of various cancer hallmarks in HNSCC. The R package “wgcna” (weighted gene co-expression network analysis) was used to construct a scale-free co-expression network. Subsequently, transcriptome profiling data and ssGSEA scores were used to screen the gene modules most related to glycolysis and hypoxia [16]. Gene significance (GS) is used to quantify the correlation between individual genes and ssGSEA scores of glycolysis and hypoxia. Module member represented the correlations between module characteristic genes and gene expression profiles. Using the p-value threshold of GS < 0.0001 and the p-value of univariate Cox regression < 0.01, 97 prognostic genes most related to glycolysis and hypoxia were identified [17]. Next, random forest was used to rank the importance of genes, and the top ten important genes were selected. The gene signature with the smaller number of genes and the more significant P value was selected from multiple combinations of ten genes and used to construct a survival model.

Survival analysis based on the risk score

The Z-score method was used to standardize the ssGSEA scores and risk scores [18]. The Kaplan–Meier method was used to perform survival analysis. The Cox proportional hazards regression model was used to evaluate the importance of each parameter to OS. Taking the median of risk scores as the cut-off value, HNSCC patients were divided into high- and low-risk groups, and the prognosis of the two groups were compared in the training set and validation set. The ROC curve was used to evaluate the prediction accuracy of the risk score. The t-ROC was used to evaluate the predictive ability (R package “survival-ROC”) [19]. Furthermore, a two-factor survival analysis combining risk score and cancer hallmarks including glycolysis and hypoxia was conducted to evaluate the impact of risk score, glycolysis and hypoxia on the prognosis of patients with HNSCC.

Comparing the survival prediction ability among various gene signatures in HNSCC

Previous studies had constructed some gene signatures that could be used to predict the survival of HNSCC patients. For example, an eight-gene signature constructed by Baoling Liu et al. and a ten-gene signature constructed by Zhaoyi Lu et al. showed good predictive ability in HNSCC [20, 21]. Therefore, in order to further evaluate the prognostic value of the genetic features constructed in this study, we compared the predictive power of the three gene signatures in HNSCC.

Drug susceptibility analysis of potential inhibitors of nine hub genes

To explore potential inhibitors available for nine hub genes, we performed drug sensitivity analysis. Drug sensitivity data for gene inhibitors were downloaded from the CellMiner™ database (version: 2020.3, database: 2.4.2, https://discover.nci.nih.gov/cellminer/home.do) [22]. The R packages “impute”, “limma”, “ggplot2” and “ggpubr” were used for data processing and visualization.

Establishment and evaluation of nomogram for survival prediction in HNSCC patients

Nomogram is an effective method for predicting the prognosis of cancer patients by simplifying complex statistical prediction models into contour maps that assess the probability of individual patients’ OS [23]. In this study, we constructed a nomogram based on the nine-gene signature to assess the probability of OS for HNSCC patients at 1-, 3-, and 5-year. Meanwhile, the predicted probability of the nomogram was compared with the observed actual probability through the calibration curve to verify the accuracy of the nomogram. Overlap with the reference line indicates that the model is accurate. In addition, t-ROC analysis was used to evaluate the survival prediction ability of this nomogram.

Correlation analysis between risk score and tumor immune microenvironment (TIME)

TIME is the immune-related complex environment in which tumor cells grow. Therefore, we conducted a correlation analysis between the risk score and TIME. Firstly, multiple analysis methods including TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL and EPIC were used to analyze the correlations between risk scores and immune cell infiltration status [24]. The immune cells involved in the analysis included CD4 T cells, CD8 T cells, B cells, neutrophils, and so on. Subsequently, the correlations between the risk score and immune cells, and the correlations between the risk score and the expression levels of immune-related molecules were explored respectively. Furthermore, the correlations between the risk score and genomic changes including gene mutation rate, gene mutation load, gene microsatellite instability and TP53 mutation were analyzed respectively.

Bioinformatics and statistical analysis

The hallmark gene sets of glycolysis and hypoxia from the MSigDB were used for GSEA analysis to verify the glycolysis and hypoxia status in the high-risk score group [25]. IBM SPSS Statistics 20 (IBM Corp., Armonk, NY, USA) and R software (version 3.5.2, http://www.r-project.org) were used to analyze data and draw figures. Z-scores were used to standardize ssGSEA scores and risk scores. The Kaplan–Meier method was used to draw the survival curve, and the log-rank test was used to assess the difference. The “wilcox.test” function was used to compare the risk scores between groups. The Cox proportional hazards regression model was used to evaluate the importance of each parameter to OS.

Results

Schematic diagram of research design

Figure 1 is a flowchart of the entire work of this research. The detailed process of constructing the OS prediction model for HNSCC patients is as follows. Firstly, glycolysis and hypoxia were identified as the main risk factors for the OS of patients with HNSCC among virous cancer hallmarks. Then, ssGSEA algorithm, WGCNA and univariate Cox regression analysis were used to screen promising candidates. Next, random forest algorithm and multiple combined screening methods were used to establish prognostic gene signature related to glycolysis and hypoxia. Finally, the prognostic value of the risk score based on the gene signature was evaluated in the training set and two independent validation sets (GSEGSE41613 and GSE42743). The patients’ information in the TCGA and GEO cohorts is shown in Table 1.

Fig. 1
figure 1

Overall flowchart of this study. LASSO, least absolute shrinkage and selection operator; HNSCC, head and neck squamous cell carcinoma; ssGSEA, single-sample gene set enrichment analysis; ROC, receiver operating characteristic; WGCNA, weighted gene co-expression network analysis

Full size image
Table 1 Clinical information of patients with HNSCC included in this study
Full size table

Glycolysis and hypoxia were identified as the main risk factors for the OS of patients with HNSCC

According to the ssGSEA scores and OS information of cancer hallmarks in the training set, the Cox coefficient of each marker was calculated and sorted. The results of univariate analysis suggested that glycolysis and hypoxia have the greatest impact on the survival of patients with HNSCC compared with other cancer hallmarks including oxidative phosphorylation, MYC_targets, angiogenesis, adipogenesis, protein secretion and DNA repair (Fig. 2A). As shown in Fig. 2B-C, the ssGSEA Z-scores of glycolysis and hypoxic were higher in the patients who died during the follow-up period compared with that of alive patients. Subsequently, using the median of Z-scores as the cut-off value, 498 HNSCC patients in the training set were divided into high- and low-risk groups. The OS rate of the high glycolytic Z-scores group was lower than that of the low glycolytic Z-scores group (HR = 1.60, P = 0.001; Fig. 2D). Meanwhile, the OS rate of the high-hypoxic Z-score group was lower than that of the low-hypoxic Z-score group (HR = 1.46, P = 0.005; Fig. 2E).

Fig. 2
figure 2

Glycolysis and hypoxia are the main risk factors for the overall survival of patients with HNSCC. A Univariate Cox regression analysis showed that glycolysis and hypoxia were the main risk factors in various cancer hallmarks. B-C The Z Score values of patients who died during the follow-up period were significantly higher than those of alive patients during the follow-up period in glycolysis and hypoxia. D-E Kaplan–Meier analysis suggested that patients with higher Z scores exhibited poorer OS in glycolysis and hypoxia. OS, overall survival; Pl3-Akt-mTOR, phosphatidylinositol-3-kinase-Akt-mammalian target of rapamycin; TGF, transforming growth factor

Full size image

Establishment of prognostic gene signature related to glycolysis and hypoxia

Using |log2FC|> 1 and FDR < 0.05, a total of 3391 DEGs in HNSCC were identified (Fig. 3A). WGCNA was conducted using transcriptome profiling data of 3391 genes and ssGSEAZ scores of the two cancer hallmarks (including glycolysis and hypoxia) in the training set. Taking β = 3 as the optimal soft threshold to ensure the scale-free co-expression of the network, the non-gray modules were generated (Fig. 3B). Then co-expression modules related to key cancer hallmarks were constructed. And the blue module was identified as the module with higher correlation with glycolysis and hypoxia (r > 0.3, P < 0.0001; Fig. 3C). Using < 0.01 as the threshold of the P value for univariate Cox regression, 97 genes from the blue module were identified as promising candidates related to the prognosis of patients with HNSCC (Fig. 3D). Next, random forest was used to rank the importance of genes, and the top 10 relatively important genes were screened out (Fig. 3E). Subsequently, the gene combination with a smaller number of genes and a more significant P value was selected from multiple combinations of 10 genes to construct a survival risk model (Fig. 3F). Finally, nine hub genes were used to construct a prognostic model for patients with HNSCC: risk score = (-0.055) * RNF144A + 0.011 * STC1 + 0.016 * P4HA1 + (-0.042) * FMNL3 + 0.002 * ANO1 + 0.002 * BASP1 + 0.004MME + (-0.106) * PLEKHG2 + 0.003 * DKK1.

Fig. 3
figure 3

Establishment of a gene signature related to glycolysis and hypoxia. A Identification of DEGs in HNSCC. Using |log2FC|> 1 and FDR < 0.05, a total of 3391 DEGs in HNSCC were identified. B 3391 DEGs were used to construct the WGCNA network to identify non-grey modules. C Construction of co-expression modules related to key cancer hallmarks. The blue module was identified as the module with higher correlation with glycolysis and hypoxia (r > 0.3, P < 0.0001). D Using univariate Cox analysis, 97 candidates related to the prognosis of patients with HNSCC were identified from the genes of the blue module (P < 0.05). E Using random forest, the top 10 genes with the highest gene importance were screened out. F The gene combination with a relatively small number of genes and a relatively significant P value was selected from the multiple combinations of 10 genes to construct a survival prediction model. DEGs, differentially expressed genes. FDR, false discovery rate. WGCNA, weighted gene co-expression network analysis

Full size image

Risk score was an independent risk factor for OS in the training set

Correlation analysis suggested the expression levels of nine hub genes were related with the Z-scores of key cancer hallmarks (including glycolysis and hypoxia) in the training set (Fig. 4A). The risk scores of the patients who died during the follow-up period were significantly higher than that of alive patients (Fig. 4B). Kaplan–Meier analysis showed that the high-risk score group exhibited worse overall survival (P < 0.001, Fig. 4C). Principal component analysis suggested that risk score could be used as a new dimension to assess the prognosis of HNSCC patients (Fig. 4D). The ROC curve showed that the AUCs of risk scores for predicting 1-, 2-, 3-, 4-, and 5-year survival rates were 0.673, 0.691, 0.717, 0.710, and 0.661, respectively, suggesting that risk score was a good model for predicting the survival of HNSCC patients (Fig. 4E). Subsequent tROC analysis showed that the survival predictive power of the risk score was better than that of other clinicopathological characters (Fig. 4F). Furthermore, univariate and multivariate Cox regression analysis showed that the risk score (HR = 1.33, p < 0.001), N stage (HR = 1.24, p < 0.01), M stage (HR = 4.25, p < 0.01) and age (HR = 1.02, p < 0.001) were independent risk factors affecting OS of HNSCC patients (Fig. 4G).

Fig. 4
figure 4

The risk score predicts poor survival in the training set. A Correlation analysis between expression levels of nine hub genes and Z-scores of key cancer hallmarks (including glycolysis and hypoxia) in the training set. B Patients who died during follow-up had a higher risk score than those who were alive. C Kaplan–Meier analysis showed that the high-risk score group exhibited worse overall survival. D Principal component analysis suggested that risk score could be used as a new dimension to evaluate the prognosis of HNSCC patients. E The ROC curve showed the prediction efficiency of the risk score for the survival rate in the training set (AUC > 0.66). F The tROC analysis showed that the predictive power of risk score was significantly higher than that of other clinical characteristics. G Univariate and multivariate Cox regression analysis showed that risk score was an independent risk factor for OS in patients with HNSCC. HR, hazard ratio; OS, overall survival; tROC, time-dependent receiver operating characteristics

Full size image

Validating the prognostic values of the nine-gene signature in independent HNSCC datasets

In order to confirm the robustness of the prognostic nine-gene signature related to glycolysis and hypoxia, two independent external cohorts (GSE41613 and GSE42743) from the GEO database were used as validation sets. The analysis results in the GSE41613 HNSCC cohort were shown in Fig. 5. Glycolysis and hypoxia pathways were significantly enriched in the high-risk group (Fig. 5A). Meanwhile, the risk score of patients who died during the follow-up period was significantly higher than that of alive patients (P < 0.0001; Fig. 5B). Kaplan–Meier analysis showed that patients with high-risk scores had worse OS than those with low-risk scores (p < 0.001; Fig. 5C). Next, the principal component analysis suggested that risk score could be used as a new dimension to assess the prognosis of HNSCC patients (Fig. 5D). The ROC curve showed that AUCs of risk scores for predicting 1-, 2-, 3-, 4-, and 5-year survival rates were 0.791, 0.778, 0.781, 0.789, and 0.766, respectively, indicating that the risk score was a good predictive model for HNSCC patients (Fig. 5E). The tROC analysis showed that the survival predictive ability of the risk score was significantly higher than other clinicopathological characteristics (Fig. 5F). In addition, multivariate Cox regression analysis suggested that risk score was an independent risk factor for PFS in HNSCC patients (HR = 2.71, p < 0.001; Fig. 5G).

Fig. 5
figure 5

Validation of the risk model in the GSE41613 dataset. A GSEA analysis showed that glycolysis and hypoxia pathways were significantly enriched in the high-risk group. B The risk score of patients who died during follow-up was higher than that of alive patients. C Kaplan–Meier analysis showed that patients with higher risk scores had worse overall survival. D Principal component analysis suggested that risk score could be used as a new dimension to evaluate the prognosis of HNSCC patients. E The ROC curve showed the prediction efficiency of the risk score in the validation dataset for the survival rate of patients with HNSCC (AUC > 0.76). F The tROC analysis showed that the survival predictive ability of risk score was significantly higher than that of other clinical features. G Univariate and multivariate Cox regression analysis showed that risk score was an independent risk factor for PFS in patients with HNSCC. HR, hazard ratio; PFS, progression-free survival; tROC, time-dependent receiver operating characteristics

Full size image

Similarly, the analysis results in the GSE42743 HNSCC cohort were shown in Fig. 6. Glycolysis and hypoxia pathways were significantly enriched in the high-risk group (Fig. 6A). The risk score of patients who died during the follow-up period was significantly higher than that of alive patients (P < 0.0001; Fig. 6B). Kaplan–Meier analysis showed that patients with high-risk scores had worse OS than those with low-risk scores (p < 0.001; Fig. 6C). Principal component analysis suggested that risk score could be used as a new dimension to assess the prognosis of HNSCC patients (Fig. 6D). The ROC curve showed that AUCs of risk scores for predicting 1-, 2-, 3-, 4-, and 5-year survival rates were 0.834, 0.920, 0.755, 0.757, and 0.853, respectively, indicating that the risk score was a good predictive model for HNSCC patients (Fig. 6E). The tROC analysis showed that the survival predictive ability of risk score was better than other clinicopathological characteristics (Fig. 6F). Moreover, multivariate Cox regression analysis showed that risk score was an independent risk factor for PFS in HNSCC patients (HR = 1.66, p < 0.001; Fig. 6G).

Fig. 6
figure 6

Validation of the risk model in the GSE42743 dataset. A GSEA analysis showed that glycolysis and hypoxia pathways were significantly enriched in the high-risk group. B The risk score of patients who died during follow-up was higher than that of alive patients. C Kaplan–Meier analysis showed that patients with higher risk scores had worse overall survival. D Principal component analysis suggested that risk score could be used as a new dimension to evaluate the prognosis of HNSCC patients. E The ROC curve showed the prediction efficiency of the risk score in the validation dataset for the survival rate of patients with HNSCC (AUC > 0.75). F The tROC analysis showed that the survival predictive ability of risk score was significantly higher than that of other clinical features. G Univariate and multivariate Cox regression analysis showed that risk score was an independent risk factor for PFS in patients with HNSCC. HR, hazard ratio; PFS, progression-free survival; tROC, time-dependent receiver operating characteristics

Full size image

Correlation between risk score and cancer hallmarks and corresponding two-factor survival analysis

As shown in Fig. 7A, the proportion of patients with high hypoxia, high glycolysis and high risk scores among patients who died during follow-up was higher, while the proportion of patients with low hypoxia, low glycolysis, and low risk scores among alive patients was higher. As shown in Fig. 7B, the Z-socres of hypoxia and glycolysis were both higher in patients with high-risk groups. Subsequently, a two-factor survival analysis combining risk score and cancer hallmarks suggested that HNSCC patients with low risk score & low glycolysis or hypoxia Z-socres showed the best OS, while HNSCC patients with high risk score & high glycolysis or hypoxia Z-socres showed the worst prognosis (Fig. 7C-D).

Fig. 7
figure 7

Two-factor survival analysis combining cancer hallmarks and risk scores. A Correlation analysis of cancer hallmarks (including hypoxia and glycolysis), risk scores and survival status of patients with HNSCC. B The Z-Scores of high-risk patients were significantly higher than those of low-risk patients both in hypoxia and glycolysis. C A two-factor survival analysis combining hypoxia and risk score suggested that high-hypoxia & high-risk score predicted a worse prognosis. D A two-factor survival analysis combining glycolysis and risk scores suggested that high-glycolysis & high-risk scores predicted a worse prognosis

Full size image

The risk score was an indicator of poor prognosis in various subgroup cohorts

As shown in Fig. 8A-H, the risk score based on the nine-gene signature could distinguish high-risk patients with poor prognosis in various subgroups divided by clinicopathological characteristics including age, gender, grades and pathological stages (p < 0.001). Meanwhile, we divided HNSCC patients into four subgroups based on TP53 and TTN mutational status, and then analyzed the effect of risk score on the survival of HNSCC patients within each subgroup. As shown in Figure S2, the risk score based on the nine-gene signature could distinguish patients with poor prognosis in four subgroups divided by TP53 and TTN gene mutations (p < 0.001).

Fig. 8
figure 8

The risk score based on the nine-gene signature is a valuable marker for poor prognosis in various subgroups divided by clinicopathological characteristics. A-H The nine-gene signature could distinguish high-risk patients in a variety of subgroups divided by clinicopathological characteristics including age, gender, grades, and pathological stages. HR, risk ratio

Full size image

This nine-gene signature was a good model for predicting the survival of HNSCC patients

Comparing the predictive power of the nine-gene signature built in this work with that of gene signatures built by other previous studies could help to further evaluate the prognostic value of the nine-gene signature. Therefore, we included two previously established gene signatures into the analysis. Included three gene signatures were as follows: the nine-gene signature was built in this work; the eight-gene signature (CBX3, GNA12, P4HA1, PLAU, PPL, RAB25, EPHX3, and HLF) was built by Baoling Liu et al.; the ten-gene signature (IL2RA, CCL5, SLC2A6, PTX3, PDGFA, INHBA, HS3ST1, TGFB1, GAS7, and RA114) was built by Zhaoyi Lu et al. As shown in Fig. 9, the predictive ability of the nine-gene signature was better than that of the other two gene signatures within eight years. Moreover, the predicted results of the nine-gene signature were basically consistent with the survival results actually observed. These results indicated that the nine-gene signature was a good prognostic model for HNSCC patients.

Fig. 9
figure 9

Comparison of survival prediction ability among various prognostic gene signatures in HNSCC. The tROC and C-index analyses were used to compare the predictive power of various gene signatures. Three gene signatures included were defined as follows: the nine-gene signature were constructed in our study; the eight-gene signature (CBX3, GNA12, P4HA1, PLAU, PPL, RAB25, EPHX3, and HLF) was constructed by Baoling Liu et al.; the ten-gene signature (IL2RA, CCL5, SLC2A6, PTX3, PDGFA, INHBA, HS3ST1, TGFBI, GAS7, and RAl14) was constructed by Zhaoyi Lu et al. A The tROC analysis showed that the predictive ability of the nine-gene signature was better than that of the other two gene signatures within eight years. B C-index analysis indicated that the prediction results of the nine-gene signature were basically consistent with the survival results actually observed in HNSCC

Full size image

Drug susceptibility analysis of related inhibitors of nine hub genes

To explore potential inhibitors of nine hub genes, we performed a drug sensitivity analysis using the CellMiner™ database. The results showed that RNF144A expression was positively correlated with the drug sensitivity of Vemurafenib, Dabrafenib, Encorafenib and ABT-199; BASP1 expression was negatively correlated with drug sensitivity of Tamoxifen, Vemurafenib, Nilotinib and Selumetinib; BASP1 expression was positively correlated with drug sensitivity of Dasatinib. STC1 expression was positively correlated with the drug sensitivity of Olaparib, Simvastatin and Rapamycin; PLEKHG2 was negatively correlated with drug sensitivity of Palbociclib, P4HA1 was negatively correlated with drug sensitivity of 6-THIOGUANINE. DKK1 was negatively correlated with drug sensitivity of Sunitinib (Figure S3) 。

Building a nomogram to predict OS in HNSCC patients

Combining the nine-gene signature, patient age and pathological stage, a comprehensive nomogram that can be used in clinical practice was constructed (Fig. 10A). Subsequently, the predictive power of this nomogram was verified by the calibration curve (Fig. 10B). Moreover, the tROC curve confirmed the good survival prediction ability of this nomogram (AUC > 0.68, Fig. 10C).

Fig. 10
figure 10

Building a nomogram based on the nine-gene signature for HNSCC patients. A A nomogram combining the nine-gene signature, M stage, N stage and age was constructed. B The calibration chart showed that the predicted 1- 3-, and 5-year survival probabilities were basically consistent with actual observations. C The tROC analysis indicated that this comprehensive nomogram had good survival prediction ability

Full size image

Correlation analysis between risk score and tumor immune infiltration

TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, EPIC and other methods were used to measure the proportion of immune cells in HNSCC. The heat map was used to explore the correlation between risk score and immune cell infiltration. The results showed that the degree of immune cell infiltration in the low-risk score group was higher than that in the high-risk score group (Fig. 11A). In addition, the expression patterns of immune cells and immune-related molecules in the high- and low-risk groups were explored. The results showed that the risk score was correlated with the immune infiltrations of B cells, CD8 + T cells, DCs, Mast cells, neutrophils, pDCs, T helper cells, Tfh, Th2 cells, TIL and Treg (Fig. 11B). Meanwhile, the risk score was significantly correlated with the expression levels of immune-related molecules such as APC co inhibition, CCR, check-point, cytolytic activity, HLA, inflammation promoting, T cell co-inhibition, T cell co-stimulation and type II IFN response (Fig. 11C). In addition, we explored the impact of TP53 and TTN mutations on immune infiltration status. As shown in Figure S4A, in the TP53 group, the mutant group exhibited a lower degree of immune infiltration compared with the wild group. Meanwhile, as shown in Figure S4B, in the TTN group, the immune infiltration status of the mutant group and the wild group was not significantly different.

Fig. 11
figure 11

Correlation analysis of risk score and tumor immune infiltration in HNSCC. A The heat map showed the correlation between the risk score and immune cell infiltration. B-C The ssGSEA was used to calculate the degrees of immune cell infiltration and the expression patterns of immune-related molecules in the high- and low-risk groups. The results suggested that the degrees of immune cell infiltration and the expression levels of immune-related molecules in the low-risk group were higher than those of the high-risk group

Full size image

Comprehensive analysis of genomic changes between high- and low-risk groups

Using genetic mutation data from HNSCC patients in TCGA, genomic changes between high- and low-risk groups were compared. The results suggested that the high-risk group had a higher rate of gene mutations (Fig. 12A-B). At the same time, the risk score of the TP53 mutant group was higher than that of the TP53 wild group (P < 0.01, Fig. 12C). And the high-risk group had a higher gene mutation load than that of the low-risk group (P < 0.05, Fig. 12D). In addition, the MSI of the high-risk group was higher than that of the low-risk group (P < 0.01, Fig. 12E).

Fig. 12
figure 12

Comprehensive comparison of genomic changes between high- and low-risk groups. A-B Comparison of gene mutation rates between the high- and low-risk groups. C Comparison of the risk scores between the TP53 mutation and the TP53 wild groups. D Comparison of the mutation load status between in the high- and low-risk groups. E Comparison of the MSI status between the high- and low-risk groups. MSI, microsatellite instability

Full size image

Discussion

HNSCC is an aggressive malignant tumor with high incidence and poor prognosis. About 60% of HNSCC patients are already in advanced stages of local tumor when they are first diagnosed. After the failure of first-line treatment, the survival of patients dropped rapidly to about 3 months. Therefore, early diagnosis and treatment are of great significance to improve the clinical efficacy and survival of HNSCC patients. The insufficient predictive power of TNM staging shows the importance of screening new prognostic biomarkers for HNSCC. Therefore, this study aims to construct prognostic gene signature based on cancer hallmarks most relevant to the prognosis of HNSCC by comprehensive bioinformatics methods. Firstly, using ssGSEA and Cox-PH regression models, glycolysis and hypoxia were identified as cancer hallmarks most related to OS in HNSCC patients. All genes of the gene modules most related to glycolysis and hypoxia were extracted for further analysis. Then, WGCNA, COX univariate regression analysis, random forest algorithm and multiple combinations were used to construct this nine-gene signature (including RNF144A, STC1, P4HA1, FMNL3, ANO1, BASP1, MME, PLEKHG2 and DKK1) with prognostic value for HNSCC patients. Next, the prognostic value of the risk score based on the nine-gene signature was validated in the training set and the validation set. In addition, the correlations between the nine-gene signature and tumor immune infiltration profiles and genome changes were explored, respectively. Therefore, this nine-gene signature is an independent prognostic predictor of patients with HNSCC.

Previous studies revealed that biomarkers including mRNA, microRNA, lncRNA, DNA methylation and protein can contribute to the early diagnosis of HNSCC, the prediction of treatment response and the early monitoring of tumor recurrence [26]. For example, SERPINE1, PLAU and ACTA1 can be used as prognostic biomarkers of HNSCC [27]. PLA, CLDN8 and CDKN2A are also prognostic biomarkers for patients with HNSCC [28]. EGFR, CDK6 or CDK4 are associated with poor prognosis in HNSCC [29]. The overexpression of TMEM16A is related to the occurrence, proliferation and migration of tumor cells, and can be used as a potential biomarker of HNSCC [30, 31]. HILPDA, CD24, TCF3, SERPINE1, INHBA, P4HA2 and ACTN1 can be used to predict the response of locally advanced HPV-negative HNSCC patients receiving postoperative chemoradiation [32]. Moreover, miR-6508-5p, miR-210-5p, miR-4306 and miR-7161-3p are independent prognostic factors for HPV-negative HNSCC [33]. And miR-99a, miR-31, miR-410, miR-424 and miR-495 help predict the radiotherapy response of HNSCC patients [34]. FAM135B methylation is also an independent prognostic biomarker of HNSCC [35]. In addition, long non-coding RNAs such as AC024592.9, LINC00941, LINC01615 and MIR9-3HG, are associated with the prognosis of HNSCC [36]. IDO1 methylation can be used to predict the response of HNSCC patients to immune checkpoint inhibitors [37]. CDKN2A methylation is involved in the occurrence, progression and metastasis of HNSCC, and it is a potential diagnostic and prognostic biomarker for patients with HNSCC [38]. In addition, increased neutrophil-to-lymphocyte ratio is significantly associated with poor OS and DSS in HNSCC [39]. Besides, a previous study reported that Wang J. et al. explored the prognostic impacts of mutated genes in HNSCC and constructed a 6-gene risk model based on the hub genes [40]. The method used in our study to construct the prognostic model was different from that used by Wang J. et al. In our study, firstly, we explored the cancer hallmarks most associated with the prognosis of HNSCC patients using transcriptome profiling data and hallmark gene sets in the Molecular Signatures Database. The results showed that glycolysis and hypoxia were screened as major risk factors for the survival of HNSCC patients. Next, we used WGCNA to identify co-expression modules related to glycolysis and hypoxia. This method of screening cancer hallmarks-related genes based on network interaction is more in line with the actual situation of network regulation in the human body. Subsequently, we constructed a prognostic model based on univariate Cox analysis, random forest and multiple combinatorial screening. Therefore, this nine-gene signature which is significantly related to glycolysis and hypoxia adds new contents to the prognostic biomarkers of HNSCC.

The prognostic nine-gene signature established in this study included RNF144A, STC1, P4HA1, FMNL3, ANO1, BASP1, MME, PLEKHG2 and DKK1. As far as we know, DKK1 is involved in GPCR signal transduction and WNT ligand antagonist negative regulation of TCF-dependent signal transduction. PLEKHG2 participates in GPCR signal transduction. MME is related to peptidase activity and endopeptidase activity. BASP1 is related to the specific binding of protein domains. ANO1 is related to the activity of calcium-activated chloride channels in cells. FMNL3 is related to malignant melanoma and is involved in pathways including GPCR signal transduction and mitotic pre-metaphase. P4HA1 is related to oxidoreductase activity. STC1 is related to diseases such as fibrosarcoma and participates in ectodermal differentiation. RNF144A participates in pathways including protein metabolism and protein ubiquitination.

Previous studies suggest that increased DKK1 expression is associated with poor prognosis in HNSCC patients [41,42,43,44]. DKK1 methylation has prognostic value in HNSCC [45]. Increased expression of ANO1 could be used as a biomarker for distant metastasis and prognosis in HNSCC [46,47,48,49]. Increased expression of P4HA1 is associated with poor prognosis and increased risk of recurrence in HNSCC patients [21, 50,51,52]. P4HA1 plays a role in promoting tumor progression [53]. As one of the glycolysis-related genes, STC1 is also associated with the prognosis of HNSCC [54]. On the other hand, the impacts of the five hub genes (including RNF144A, FMNL3, BASP1, MME and PLEKHG2) on the prognosis of HNSCC patients have not yet been reported. Therefore, these results suggest that the mechanism of action of these nine hub genes in HNSCC deserves to be further explored.

In this study, we constructed a new nine-gene signature related to cancer hallmarks (including glycolysis and hypoxia). Survival analysis results confirmed that this nine-gene signature including RNF144A, STC1, P4HA1, FMNL3, ANO1, BASP1, MME, PLEKHG2 and DKK1 had good prognostic value for HNSCC patients. It is worth mentioning that recent studies had established some prognostic gene signatures for HNSCC patients as a supplement to traditional pathological staging. Comparing the survival prediction ability of previously reported gene signatures with the nine-gene signature will contribute to further assessing the prognostic values of these gene signatures. Therefore, we compared the survival prediction ability between the nine-gene signature constructed in this study and the previously reported two gene signatures. The results showed that the predictive ability of the nine-gene signature was better than the other two gene signatures within eight years, suggesting that this nine-gene signature was a good prognostic model in HNSCC.

The results of this study suggest that the risk score based on nine-gene signature is related to the immune cells infiltration status, the expression profiles of immune-related molecules and genome changes in HNSCC. Previous studies have shown that the immune microenvironment of HNSCC is characterized by changes in immune cells and immune checkpoints that make the balance of the immune environment beneficial to immune suppression, thereby allowing tumor cells to escape immune surveillance. Therefore, immunotherapy is expected to be a potentially beneficial supplement to the standard treatment of HNSCC [55,56,57,58]. Recent studies have shown that the immune infiltrating status is with prognostic value in HNSCC [59]. The prognosis of HPV-positive HNSCC patients is better than that of HPV-negative HNSCC patients, possibly due in part to the enhanced immune activation of the tumor microenvironment in HPV-positive HNSCC [60]. In addition, programmed death ligand 1 (PD-L1) is an immune checkpoint mainly located on the surface of tumor cells, and the positive expression of PD-L1 is associated with a better prognosis of HNSCC patients [61, 62]. Therefore, immunotherapies, especially those targeting the PD1 receptor or its ligand PD-L1, have shown significant efficacy in HNSCC [63]. In terms of genomic changes, tumor mutation burden (TMB) has been considered as a predictor of immune checkpoint inhibitors (ICIs) response. High TMB can identify HNSCC patients with poor prognosis after concurrent radiotherapy and chemotherapy [64]. PD-1 overexpression is associated with a good prognosis [65]. Currently known biomarkers predicting response to immune checkpoint inhibitors include PD-L1 expression, human papilloma virus infection, and microsatellite instability [66]. Therefore, the nine-gene signature constructed in this study were related to the degree of immune cell infiltration, the expression levels of immune-related molecules, and genomic changes. Meanwhile, we believe that two aspects should be considered to understand the results of the immune infiltration analysis. On the one hand, the predominant immune cells exhibited a higher degree of immune infiltration in the low-risk group. This may make the low-risk group have better prognosis. On the other hand, as shown in Fig. 11B, some immunosuppressive cells (such as T-reg cells) also had higher expression in the low-risk group. In addition, as shown in Fig. 11C, some immune-related molecules with opposite functions (such as T-cell co-stimulatory and T-cell co-inhibitory) were over-expressed at the same time in the low-risk group. We think that this is due to the complexity of the microenvironment. For example, Kubli, SP et al. proposed that the complexity of immune cell-cancer cell interaction is an important reason for the failure of tumor immunotherapy [67]. Besides, Feng, C et al., Prokhnevska, N et al. and Wierz, M et al. described the complexity of the tumor microenvironment in bladder cancer, prostate cancer, and chronic lymphocytic leukemia, respectively [68,69,70]. Taken together, we believe that the results of the immune infiltration analysis not only showed that the main immune cells exhibited a higher degree of immune infiltration in the low-risk group, but also provide some clues about the complexity of immune cell roles in HNSCC. These are the reasons why the immune infiltration status in HNSCC deserves to be further explored.

There are some limitations in this study. This study used a method of interactive verification among multiple independent datasets to verify the prognostic significance of the nine-gene signature. However, experimental validation is still an important step to further verify the prognostic value of this model. And the lack of experimental verification is the limitation of this study. In addition, this is a retrospective study, so the robustness and clinical usefulness of this nine-gene signature needs to be verified in prospective clinical trials.

In conclusion, this study identified a new prognostic gene signature related to glycolysis and hypoxia, including RNF144A, STC1, P4HA1, FMNL3, ANO1, BASP1, MME, PLEKHG2 and DKK1. Besides, a nomogram based on the nine-gene signature was constructed for clinical practice. This nine-gene signature can not only be used as a prognostic biomarker to help clinicians develop more personalized treatments for HNSCC patients, but also is expected to become a potential therapeutic target for HNSCC.

Availability of data and materials

The datasets analyzed in the current study are available in the TCGA repository (http://cancergenome.nih.gov/) and the GEO (https://www.ncbinlm.nih.gov/geo/).

References

  1. Cohen N, Fedewa S, Chen AY. Epidemiology and Demographics of the Head and Neck Cancer Population. Oral Maxillofac Surg Clin N Am. 2018;30(4):381.

    Article  Google Scholar 

  2. Li H, Torabi SJ, Yarbrough WG, Mehra S, Osborn HA, Judson B. Association of Human Papillomavirus Status at Head and Neck Carcinoma Subsites With Overall Survival. JAMA Otolaryngol-Head Neck Surg. 2018;144(6):519–25.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Budach V, Tinhofer I. Novel prognostic clinical factors and biomarkers for outcome prediction in head and neck cancer: a systematic review. Lancet Oncol. 2019;20(6):E313–26.

    Article  PubMed  Google Scholar 

  4. Kumar D, New J, Vishwakarma V, Joshi R, Enders J, Lin FC, Dasari S, Gutierrez WR, Leef G, Ponnurangam S, et al. Cancer-Associated Fibroblasts Drive Glycolysis in a Targetable Signaling Loop Implicated in Head and Neck Squamous Cell Carcinoma Progression. Cancer Res. 2018;78(14):3769–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Fleming JC, Woo J, Moutasim K, Mellone M, Frampton SJ, Mead A, Ahmed W, Wood O, Robinson H, Ward M, et al. HPV, tumour metabolism and novel target identification in head and neck squamous cell carcinoma. Br J Cancer. 2019;120(3):356–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Wang LH, Bai J, Duan P. Prognostic value of F-18-FDG PET/CT functional parameters in patients with head and neck cancer: a meta-analysis. Nucl Med Commun. 2019;40(4):361–9.

    Article  CAS  PubMed  Google Scholar 

  7. Chen TY, Hsieh YT, Huang JM, Liu CJ, Chuang LT, Huang PC, Kuo TY, Chia HY, Chou CY, Chang CW, et al. Determination of Pyruvate Metabolic Fates Modulates Head and Neck Tumorigenesis. Neoplasia. 2019;21(7):641–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Wiedenmann N, Grosu AL, Buchert M, Rischke HC, Ruf J, Bielak L, Majerus L, Ruhle A, Bamberg F, Baltas D, et al. The utility of multiparametric MRI to characterize hypoxic tumor subvolumes in comparison to FMISO PET/CT. Consequences for diagnosis and chemoradiation treatment planning in head and neck cancer. Radiother Oncol. 2020;150:128–35.

    Article  CAS  PubMed  Google Scholar 

  9. Nicolay NH, Ruhle A, Wiedenmann N, Niedermann G, Mix M, Weber WA, Baltas D, Werner M, Kayser G, Grosu AL. Lymphocyte Infiltration Determines the Hypoxia-Dependent Response to Definitive Chemoradiation in Head-and-Neck Cancer: Results from a Prospective Imaging Trial. J Nucl Med. 2021;62(4):471–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Yan XW, Cao NQ, Chen Y, Lan HY, Cha JH, Yang WH, Yang MH. MT4-MMP promotes invadopodia formation and cell motility in FaDu head and neck cancer cells. Biochem Biophys Res Commun. 2020;522(4):1009–14.

    Article  CAS  PubMed  Google Scholar 

  11. Linge A, Schmidt S, Lohaus F, Krenn C, Bandurska-Luque A, Platzek I, von Neubeck C, Appold S, Nowak A, Gudziol V, et al. Independent validation of tumour volume, cancer stem cell markers and hypoxia-associated gene expressions for HNSCC after primary radiochemotherapy. Clin Transl Radiat Oncol. 2019;16:40–7.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Stieb S, Eleftheriou A, Warnock G, Guckenberger M, Riesterer O. Longitudinal PET imaging of tumor hypoxia during the course of radiotherapy. Eur J Nucl Med Mol Imaging. 2018;45(12):2201–17.

    Article  PubMed  Google Scholar 

  13. Wang Z, Wang Y, Yang T, Xing H, Wang Y, Gao L, Guo X, Xing B, Wang Y, Ma W. Machine learning revealed stemness features and a novel stemness-based classification with appealing implications in discriminating the prognosis, immunotherapy and temozolomide responses of 906 glioblastoma patients. Brief Bioinform. 2021;22(5):bbab032.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462(7269):108–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27(12):1739–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Li W, Lu J, Ma Z, Zhao J, Liu J. An Integrated Model Based on a Six-Gene Signature Predicts Overall Survival in Patients With Hepatocellular Carcinoma. Front Genet. 2019;10:1323.

    Article  CAS  PubMed  Google Scholar 

  18. Yang X, Weng X, Yang Y, Zhang M, Xiu Y, Peng W, Liao X, Xu M, Sun Y, Liu X. A combined hypoxia and immune gene signature for predicting survival and risk stratification in triple-negative breast cancer. Aging (Albany NY). 2021;13(15):19486–509.

    Article  CAS  Google Scholar 

  19. Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.

    Article  CAS  PubMed  Google Scholar 

  20. Lu ZY, Deng XL, Li H. Prognostic Value of a Ten-Gene Signature in HNSCC Patients Based on Tumor-Associated Macrophages Expression Profiling. Front Oncol. 2020;10:12.

    Article  Google Scholar 

  21. Liu BL, Su QP, Ma JH, Chen C, Wang LJ, Che FY, Heng XY. Prognostic Value of Eight-Gene Signature in Head and Neck Squamous Carcinoma. Front Oncol. 2021;11:13.

    Google Scholar 

  22. Reinhold WC, Sunshine M, Liu H, Varma S, Kohn KW, Morris J, Doroshow J, Pommier Y. Cell Miner: a web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the NCI-60 cell line set. Cancer Res. 2012;72(14):3499–511.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Park SY. Nomogram: An analogue tool to deliver digital knowledge. J Thorac Cardiovasc Surg. 2018;155(4):1793.

    Article  PubMed  Google Scholar 

  24. Wang M, Zhu J, Zhao F, Xiao J. Transcriptome Analyses Identify a Metabolic Gene Signature Indicative of Antitumor Immunosuppression of EGFR Wild Type Lung Cancers With Low PD-L1 Expression. Front Oncol. 2021;11:643503.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Li J, Du J, Wang Y, Jia H. A Coagulation-Related Gene-Based Prognostic Model for Invasive Ductal Carcinoma. Front Genet. 2021;12:722992.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Arantes L, De Carvalho AC, Melendez ME, Carvalho AL. Serum, plasma and saliva biomarkers for head and neck cancer. Expert Rev Mol Diagn. 2018;18(1):85–112.

    Article  CAS  PubMed  Google Scholar 

  27. Yang K, Zhang SZ, Zhang DS, Tao Q, Zhang TQ, Liu GJ, Liu XG, Zhao TD. Identification of SERPINE1, PLAU and ACTA1 as biomarkers of head and neck squamous cell carcinoma based on integrated bioinformatics analysis. Int J Clin Oncol. 2019;24(9):1030–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Zhao XY, Sun SY, Zeng XQ, Cui L. Expression profiles analysis identifies a novel three-mRNA signature to predict overall survival in oral squamous cell carcinoma. Am J Cancer Res. 2018;8(3):450–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Sayans MP, Petronacci CMC, Pouso AIL, Iruegas EP, Carrion AB, Penaranda JMS, Garcia AG. Comprehensive Genomic Review of TCGA Head and Neck Squamous Cell Carcinomas (HNSCC). J Clin Med. 2019;8(11):18.

    Google Scholar 

  30. Godse NR, Khan N, Yochum ZA, Gomez-Casal R, Kemp C, Shiwarski DJ, Seethala RS, Kulich S, Seshadri M, Burns TF, et al. TMEM16A/ANO1 Inhibits Apoptosis Via Downregulation of Bim Expression. Clin Cancer Res. 2017;23(23):7324–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Ji QS, Guo S, Wang XZ, Pang CL, Zhan Y, Chen YF, An HL. Recent advances in TMEM16A: Structure, function, and disease. J Cell Physiol. 2019;234(6):7856–73.

    Article  CAS  PubMed  Google Scholar 

  32. Schmidt S, Linge A, Zwanenburg A, Leger S, Lohaus F, Krenn C, Appold S, Gudziol V, Nowak A, von Neubeck C, et al. Development and Validation of a Gene Signature for Patients with Head and Neck Carcinomas Treated by Postoperative Radio(chemo)therapy. Clin Cancer Res. 2018;24(6):1364–74.

    Article  CAS  PubMed  Google Scholar 

  33. Hess J, Unger K, Maihoefer C, Schuttrumpf L, Wintergerst L, Heider T, Weber P, Marschner S, Braselmann H, Samaga D, et al. A Five-MicroRNA Signature Predicts Survival and Disease Control of Patients with Head and Neck Cancer Negative for HPV Infection. Clin Cancer Res. 2019;25(5):1505–16.

    Article  PubMed  Google Scholar 

  34. Chen L, Wen YH, Zhang JW, Sun W, Lui VWY, Wei Y, Chen FH, Wen WP. Prediction of radiotherapy response with a 5-microRNA signature-based nomogram in head and neck squamous cell carcinoma. Cancer Med. 2018;7(3):726–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zhou CC, Ye M, Ni SM, Li Q, Ye D, Li JY, Shen ZS, Deng HX. DNA methylation biomarkers for head and neck squamous cell carcinoma. Epigenetics. 2018;13(4):398–409.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Hu Y, Guo GY, Li JJ, Chen J, Tan PQ. Screening key lncRNAs with diagnostic and prognostic value for head and neck squamous cell carcinoma based on machine learning and mRNA-lncRNA co-expression network analysis. Cancer Biomark. 2020;27(2):195–206.

    Article  CAS  PubMed  Google Scholar 

  37. Sailer V, Sailer U, Bawden EG, Zarbl R, Wiek C, Vogt TJ, Dietrich J, Loick S, Grunwald I, Toma M, et al. DNA methylation of indoleamine 2,3-dioxygenase 1 (IDO1) in head and neck squamous cell carcinomas correlates with IDO1 expression, HPV status, patients’ survival, immune cell infiltrates, mutational load, and interferon gamma signature. EBioMedicine. 2019;48:341–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zhou CC, Shen ZS, Ye D, Li Q, Deng HX, Liu HG, Li JY. The Association and Clinical Significance of CDKN2A Promoter Methylation in Head and Neck Squamous Cell Carcinoma: a Meta-Analysis. Cell Physiol Biochem. 2018;50(3):868–82.

    Article  CAS  PubMed  Google Scholar 

  39. Tham T, Bardash Y, Herman SW, Costantino PD. Neutrophil-to-lymphocyte ratio as a prognostic indicator in head and neck cancer: A systematic review and meta-analysis. Head Neck-J Sci Spec Head Neck. 2018;40(11):2546–57.

    Article  Google Scholar 

  40. Wang JC, Chen X, Tian YX, Zhu GC, Qin YX, Chen X, Pi LM, Wei M, Liu GC, Li ZX, et al. Six-gene signature for predicting survival in patients with head and neck squamous cell carcinoma. Aging-US. 2020;12(1):767–83.

    Article  CAS  Google Scholar 

  41. Gao HH, Li LS, Xiao M, Guo YW, Shen Y, Cheng LX, Tang M. Elevated DKK1 expression is an independent unfavorable prognostic indicator of survival in head and neck squamous cell carcinoma. Cancer Manag Res. 2018;10:5083–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Hu Y, Liu MY, Xu SW, Li SK, Yang MF, Su T, Yuan ZH, Peng HW. The Clinical Significance of Dickkopf Wnt Signaling Pathway Inhibitor Gene Family in Head and Neck Squamous Cell Carcinoma. Med Sci Monitor. 2020;26:13.

    Article  Google Scholar 

  43. Liu XY, Chen JR, Lu W, Zeng ZH, Li JL, Jiang XP, Gao YP, Gong Y, Wu QJ, Xie CH. Systematic Profiling of Immune Risk Model to Predict Survival and Immunotherapy Response in Head and Neck Squamous Cell Carcinoma. Front Genet. 2020;11:15.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Zhang Y, Chen P, Zhou Q, Wang HY, Hua QQ, Wang J, Zhong HL. A Novel Immune-Related Prognostic Signature in Head and Neck Squamous Cell Carcinoma. Front Genet. 2021;12:15.

    CAS  Google Scholar 

  45. Paluszczak J, Hemmerling D, Kostrzewska-Poczekaj M, Jarmuz-Szymczak M, Grenman R, Wierzbicka M, Baer-Dubowska W. Frequent hypermethylation of WNT pathway genes in laryngeal squamous cell carcinomas. J Oral Pathol Med. 2014;43(9):652–7.

    Article  CAS  PubMed  Google Scholar 

  46. Duvvuri U, Shiwarski DJ, Xiao D, Bertrand C, Huang X, Edinger RS, Rock JR, Harfe BD, Henson BJ, Kunzelmann K, et al. TMEM16A Induces MAPK and Contributes Directly to Tumorigenesis and Cancer Progression. Cancer Res. 2012;72(13):3270–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Ruiz C, Martins JR, Rudin F, Schneider S, Dietsche T, Fischer CA, Tornillo L, Terracciano LM, Schreiber R, Bubendorf L, et al. Enhanced Expression of ANO1 in Head and Neck Squamous Cell Carcinoma Causes Cell Migration and Correlates with Poor Prognosis. PLoS ONE. 2012;7(8):12.

    Article  Google Scholar 

  48. Reddy RB, Bhat AR, James BL, Govindan SV, Mathew R, Ravindra DR, Hedne N, Illiayaraja J, Kekatpure V, Khora SS, et al. Meta-Analyses of Microarray Datasets Identifies ANO1 and FADD as Prognostic Markers of Head and Neck Cancer. PLoS ONE. 2016;11(1):21.

    Article  Google Scholar 

  49. Lek SM, Li K, Tan QX, Shannon NB, Ng WH, Hendrikson J, Tan JWS, Lim HJ, Chen YD, Koh KKN, et al. Pairing a prognostic target with potential therapeutic strategy for head and neck cancer. Oral Oncol. 2020;111:8.

    Article  Google Scholar 

  50. Kappler M, Kotrba J, Kaune T, Bache M, Rot S, Bethmann D, Wichmann H, Guttler A, Bilkenroth U, Horter S, et al. P4HA1: A single-gene surrogate of hypoxia signatures in oral squamous cell carcinoma patients. Clin Transl Radiat Oncol. 2017;5:6–11.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Zhao CG, Zhou YR, Ma HW, Wang JH, Guo HL, Liu H. A four-hypoxia-genes-based prognostic signature for oral squamous cell carcinoma. BMC Oral Health. 2021;21(1):13.

    Article  Google Scholar 

  52. Wu ZH, Zhong Y, Zhou T, Xiao HJ. miRNA biomarkers for predicting overall survival outcomes for head and neck squamous cell carcinoma. Genomics. 2021;113(1):135–41.

    Article  CAS  PubMed  Google Scholar 

  53. Li Q, Shen ZS, Wu ZH, Shen Y, Deng HX, Zhou CC, Liu HG. High P4HA1 expression is an independent prognostic factor for poor overall survival and recurrent-free survival in head and neck squamous cell carcinoma. J Clin Lab Anal. 2020;34(3):10.

    Article  Google Scholar 

  54. Liu YC, Yin SH. A Novel Prognostic Index Based on the Analysis of Glycolysis-Related Genes in Head and Neck Squamous Cell Carcinomas. J Oncol. 2020;2020:13.

    Article  Google Scholar 

  55. Whiteside TL. Head and Neck Carcinoma Immunotherapy: Facts and Hopes. Clin Cancer Res. 2018;24(1):6–13.

    Article  CAS  PubMed  Google Scholar 

  56. Leemans CR, Snijders PJF, Brakenhoff RH. The molecular landscape of head and neck cancer. Nat Rev Cancer. 2018;18(5):269–82.

    Article  CAS  PubMed  Google Scholar 

  57. Ferris RL, Blumenschein G, Fayette J, Guigay J, Colevas AD, Licitra L, Harrington KJ, Kasper S, Vokes EE, Even C, et al. Nivolumab vs investigator’s choice in recurrent or metastatic squamous cell carcinoma of the head and neck: 2-year long-term survival update of CheckMate 141 with analyses by tumor PD-L1 expression. Oral Oncol. 2018;81:45–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Mehra R, Seiwert TY, Gupta S, Weiss J, Gluck I, Eder JP, Burtness B, Tahara M, Keam B, Kang H, et al. Efficacy and safety of pembrolizumab in recurrent/metastatic head and neck squamous cell carcinoma: pooled analyses after long-term follow-up in KEYNOTE-012. Br J Cancer. 2018;119(2):153–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Zhang XM, Song LJ, Shen J, Yue H, Han YQ, Yang CL, Liu SY, Deng JW, Jiang Y, Fu GH, et al. Prognostic and predictive values of immune infiltrate in patients with head and neck squamous cell carcinoma. Hum Pathol. 2018;82:104–12.

    Article  CAS  PubMed  Google Scholar 

  60. Chen XH, Yan BQ, Lou HH, Shen ZJ, Tong FJ, Zhai AX, Wei LL, Zhang FM. Immunological network analysis in HPV associated head and neck squamous cancer and implications for disease prognosis. Mol Immunol. 2018;96:28–36.

    Article  CAS  PubMed  Google Scholar 

  61. Yang WF, Wong MCM, Thomson PJ, Li KY, Su YX. The prognostic role of PD-L1 expression for survival in head and neck squamous cell carcinoma: A systematic review and meta-analysis. Oral Oncol. 2018;86:81–90.

    Article  CAS  PubMed  Google Scholar 

  62. Colevas AD, Bahleda R, Braiteh F, Balmanoukian A, Brana I, Chau NG, Sarkar I, Molinero L, Grossman W, Kabbinavar F, et al. Safety and clinical activity of atezolizumab in head and neck cancer: results from a phase I trial. Ann Oncol. 2018;29(11):2247–53.

    Article  CAS  PubMed  Google Scholar 

  63. Solomon B, Young RJ, Rischin D. Head and neck squamous cell carcinoma: Genomics and emerging biomarkers for immunomodulatory cancer treatments. Semin Cancer Biol. 2018;52:228–40.

    Article  CAS  PubMed  Google Scholar 

  64. Eder T, Hess AK, Konschak R, Stromberger C, Johrens K, Fleischer V, Hummel M, Balermpas P, von der Grun J, Linge A, et al. Interference of tumour mutational burden with outcome of patients with head and neck cancer treated with definitive chemoradiation: a multicentre retrospective study of the German Cancer Consortium Radiation Oncology Group. Eur J Cancer. 2019;116:67–76.

    Article  CAS  PubMed  Google Scholar 

  65. Lecerf C, Kamal M, Vacher S, Chemlali W, Schnitzler A, Morel C, Dubot C, Jeannot E, Meseure D, Klijanienko J, et al. Immune gene expression in head and neck squamous cell carcinoma patients. Eur J Cancer. 2019;121:210–23.

    Article  CAS  PubMed  Google Scholar 

  66. Saleh K, Eid R, Haddad FGH, Khalife-Saleh N, Kourie HR. New developments in the management of head and neck cancer - impact of pembrolizumab. Therap Clin Risk Manag. 2018;14:295–303.

    Article  CAS  Google Scholar 

  67. Kubli SP, Berger T, Araujo DV, Siu LL, Mak TW. Beyond immune checkpoint blockade: emerging immunological strategies. Nat Rev Drug Discov. 2021;20(12):899–919.

    Article  CAS  PubMed  Google Scholar 

  68. Feng C, Wang X, Tao Y, Xie Y, Lai Z, Li Z, Hu J, Tang S, Pan L, He L, et al. Single-Cell Proteomic Analysis Dissects the Complexity of Tumor Microenvironment in Muscle Invasive Bladder Cancer. Cancers (Basel). 2021;13(21):5440.

    Article  CAS  Google Scholar 

  69. Prokhnevska N, Emerson DA, Kissick HT, Redmond WL. Immunological Complexity of the Prostate Cancer Microenvironment Influences the Response to Immunotherapy. Adv Exp Med Biol. 2019;1210:121–47.

    Article  CAS  PubMed  Google Scholar 

  70. Wierz M, Janji B, Berchem G, Moussay E, Paggetti J. High-dimensional mass cytometry analysis revealed microenvironment complexity in chronic lymphocytic leukemia. OncoImmunology. 2018;7(8):e1465167.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Disclosures

A preprint has previously been published in “Research Square” [Liu Jun et al. 22 November 2021]: (https://doi.org/10.21203/rs.3.rs-1067602/v1 [61].

Funding

This work is not funded.

Author information

Author notes

  1. Jun Liu and Jianjun Lu are co-first authors..

Authors and Affiliations

  1. Reproductive Medicine Center, Yue Bei People’s Hospital, Shantou University Medical College, 133 Huimin South Road, Shaoguan, 512025, China

    Jun Liu & Wenli Li

  2. Medical Research Center, Yue Bei People’s Hospital, Shantou University Medical College, Shaoguan, 512025, China

    Jun Liu

  3. Department of Medical Affairs, First Affiliated Hospital of Sun Yat-Sen University, Guangzhou, 510080, China

    Jianjun Lu

  4. The Second School of Clinical Medicine, Southern Medical University, Guangzhou, 510080, China

    Jianjun Lu

Authors
  1. Jun Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jianjun Lu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Wenli Li
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Wenli Li designed the study and revised the manuscript. Jun Liu collected and analyzed the data. Jianjun Lu interpreted the data and drafted the manuscript. All authors have read and approved the final manuscript. Jun Liu and Jianjun Lu contributed equally to this work and are co-first authors.

Corresponding author

Correspondence to Wenli Li.

Ethics declarations

Ethics approval and consent to participate

This study did not involve human participants or human material. The transcriptome data and clinical information were download from the TCGA and GEO databases which were publicly available. Therefore, this study did not require the approval of the local ethics committee. All methods were performed in accordance with the declaration of Helsinki and relevant and regulations.

Consent for publication

Not applicable.

Competing interests

All authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, J., Lu, J. & Li, W. Transcriptome analysis reveals the prognostic and immune infiltration characteristics of glycolysis and hypoxia in head and neck squamous cell carcinoma. BMC Cancer 22, 352 (2022). https://doi.org/10.1186/s12885-022-09449-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-022-09449-9

Keywords

BMC Cancer

ISSN: 1471-2407

Contact us