Skip to main content
Download PDF

A prognostic Risk Score model for oral squamous cell carcinoma constructed by 6 glycolysis-immune-related genes

BMC Oral Health volume 22, Article number: 324 (2022) Cite this article

Abstract

Background

Oral squamous cell carcinoma (OSCC) is the most frequent tumor of the head and neck. The glycolysis-related genes and immune-related genes have been proven prognostic values in various cancers. Our study aimed to test the prognostic value of glycolysis-immune-related genes in OSCC.

Methods

Data of OSCC patients were obtained from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Enrichment analysis was applied to the glycolysis- and immune-related genes screened by differential expression analysis. Univariate Cox and LASSO Cox analyses were used to filtrate the genes related to the prognosis of OSCC and to construct Risk Score model.

Results

A Risk Score model was constructed by six glycolysis-immune-related genes (including ALDOC, VEGFA, HRG, PADI3, IGSF11 and MIPOL1). High risk OSCC patients (Risk Score >−0.3075) had significantly worse overall survival than that of low risk patients (Risk Score <−0.3075).

Conclusions

The Risk Score model constructed basing on 6 glycolysis-immune-related genes was reliable in stratifying OSCC patients with different prognosis.

Peer Review reports

Background

Oral squamous cell carcinoma (OSCC) is the most commonly occurred malignancy in the oral cavity [1]. Most OSCCs are correlated with oral precursor lesions [2], and the 5-year overall survival of OSCC patients is less than 40% [3]. Many factors may trigger OSCC, including tobacco smoking, excessive alcohol drinking, chewing betel quid and human papillomavirus (HPV) infection [4]. Moreover, as the most frequent subtype of head and neck squamous cell carcinoma (HNSCC), the metastasis rate of OSCC is high, besides resistance to traditional chemotherapy is usually observed in OSCC patients, leading to undesirable prognosis [5]. Novel biomarkers correlated with the prognosis of OSCC may help to stratify patients with different prognosis and design individual-specific therapy, for example FKBP51 [6] is promising to predict the prognosis of OSCC patients. Whereas, the exploration of more effective and accurate prognostic biomarkers remains a continuous challenge for the medical industry.

Glycolysis, also called Embden-Meyerhof pathway, is an essential metabolic pathway and provides anaerobic energy for body function [7]. In aerobic conditions, pyruvate from glycolysis produces adenosine triphosphate (ATP) for cellular process through oxidative phosphorylation; and in anaerobic conditions, the pyruvate undergoes anaerobic glycolysis [8]. In cells that are not able to generate adequate ATP for body function via oxidative phosphorylation, anaerobic glycolysis may be a way to produce energy [9]. But cancer cells mainly depend on the aerobic glycolysis in order to rapidly provide energy for the tumors even in the presence of sufficient oxygen, which is known as Warburg effect [10]. Therefore, tumor cells intake more glucose to perform aerobic glycolysis. The increased glucose level and the overexpression of glucose transporter proteins (HIF-1α and GLUT-1) are connected with poor prognosis of OSCC [11]. Gong et al. have recently demonstrated that PER1 is a suppressor of glycolysis in OSCC, involving the regulation of cell proliferation [12]. Several aerobic glycolysis related genes have been reported with prognostic or diagnostic values not only in OSCC [13] but also in some other types of cancers, such as breast cancer [14]. Change of ATP supply from oxidative phosphorylation to aerobic glycolysis is thought be the biomarker of T cell activation [15]. The aerobic glycolysis promotes activation of T cells via phosphoinositide 3-kinase (PI3K)/Akt signaling [16]. The activated T cells produce more lactate via increasing lactate dehydrogenase A (LDHA) to support aerobic glycolysis, and increased LDHA is also implicated in poor prognosis of cancer patients [17]. Aerobic glycolysis, a well-known resistance factor for anticancer therapies, is associated with the immunotherapy for cancer patients. Aerobic glycolysis can impact the tumor immunosuppression via a network of pathways in breast cancer [18]. The prognosis of cancers are also closely related to the immune microenvironment and novel prognostic biomarkers for example MIR155HG [19] are reported to have correlation with immune infiltration. However, the potential impacts of glycolysis and immune have been seldom studied in OSCC.

Accordingly, basing on the publicly obtained data, the focus of our work is to discuss the prognostic value of glycolysis-immune-related genes for OSCC patients and to construct a prognostic model using glycolysis-immune-related genes for separating OSCC patients with different prognosis.

Methods

Data sources

We downloaded data of 306 OSCC patients, including mRNA expression profiles and corresponding clinical information, and listed the clinical information in Table 1. The data were collected from the Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/). Moreover, we downloaded datasets comprised of gene expression profiles and complete survival information of 246 OSCC patients from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database and numbered them GSE85446 (66), GSE65858 (83) and GSE41613 (97). The data of OSCC patients were obtained by the Agilent-014850 Whole Human Genome Microarray 4×44 K G4112F, Illumina HumanHT-12 V4.0 expression beadchip and Affymetrix Human Genome U133 Plus 2.0 Array. Datasets GSE85446, GSE65858, and GSE41613 were merged as meta-GEO dataset, using for subsequent validation.

Table 1 Clinicopathological characteristics of OSCC patients from TCGA database
Full size table

Differential expression analysis

Limma package [20] in R programming software (version 4.1.0, the same below) was applied to the data we collected above to identify differentially expressed genes (DEGs). The |Log2FC|> 0.7 and multiple testing adjusted p value < 0.05 were used as threshold.

Functional enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) enrichment analyses were applied to the DEGs using “clusterProfiler” package [21] in R programming software. The GO terms and KEGG pathways at the significant level (p value < 0.05, adjusted by Benjamini and Hochberg method) were employed.

Cluster analysis and calculation of ImmuneScore

The mRNA expression profiles were subjected to cluster analysis based on the “K-mean” method utilizing R programming software. The ImmuneScore of each patient was calculated by the “ESTIMATE” function.

LASSO Cox regression analysis

Univariate Cox analysis was applied to expression profiles of the OSCC patients based on the expression values of selected genes. To screen the genes which were significantly associated with the overall survival (OS) of OSCC, p value < 0.01 was used as threshold. To further optimize the genes significantly associated with OS of OSCC, LASSO Cox regression analysis was performed on the screen DEGs using glmnet package [22] in R programming software. Risk Score of each patient was calculated by the following formula:

$${\text{Risk Score}} = \sum_{{{\text{i}} = {1}}}^{{\text{n}}} {\text{Coef}}_{{\text{i}}} {\text{*X}}_{{\text{i}}} {,}$$

Coefi was the LASSO Cox risk coefficient and Xi was the expression value of genes (mRNA expression in this research). Risk Score was tested using the survival, survminer and two-sided log-rank test in R programming software. Then patients were assigned into high-risk and low-risk groups according to the median of Risk Score.

Kaplan–Meier survival analysis

OS rates of OSCC patients were estimated using the survival and survminer packages in R programming software. The significance of difference of OS rates between high- and low-risk groups was tested by the log-rank or breslow test. Multivariate Cox regression model was used to analyze the independence of the prognostic value of Risk Score.

Proportion of immune cell infiltration

We used “CIBERSORT” [23] software to estimate the infiltrations of immune cells for OSCC patients. The “CIBERSORT” software, utilizing the deconvolution algorithm, was based on the gene expression matrix. And it characterized the composition of the immune cells by the predisposed 547 barcode genes. The sum of the estimated proportion of the 22 immune cells was one. The significance of difference of immune infiltration ratios was tested by the wilcoxn method.

Construction of nomogram model

We used rms package in R programming software to construct a prognostic nomogram model by using all independent prognostic factors examined by multivariate Cox regression model and examined the efficiency of the nomogram model by drawing a calibration curve of nomogram.

Results

Patients with different prognosis defined by glycolysis-immune-related genes

Cluster analysis was applied to the patients’ mRNA expression profiles utilizing the 296 glycolysis-related genes downloaded from the GSEA database (Additional file 1: Table S1). According to the sum of the square errors (SSE) (Fig. 1A), the patients were divided into 2 clusters (k = 2) (Fig. 1B). The OS between cluster1 and cluster2 differed significantly via Kaplan–Meier survival analysis (Fig. 1C). The prognosis of cluster1 (low-glycolysis group) was significantly better than cluster2 (high-glycolysis group) (p = 0.018).

Fig. 1
figure 1

Different prognosis defined by glycolysis in combination with immunity. A The elbow figure to determine the number of clusters. The x-axis and y-axis are the number of cluster K and the sum of the square errors. K = 2 is determined by the slope of the curve. B Clustering sketch of patients. Different colors represent different clusters. C The Kaplan–Meier survival curves defined by glycolysis. The x-axis and y-axis are time and survival probability. The blue is cluster1 (low-glycolysis), and the red is cluster 2 (high-glycolysis). The p value is calculated by the log-rank test. D The best cutoff value of the two groups defined by ImmuneScore. E The Kaplan–Meier survival curves defined by ImmuneScore. The x-axis and y-axis are time and survival probability. The blue is high-ImmuneScore, and the red is low-ImmuneScore. F The Kaplan–Meier survival curves defined by glycolysis in combination with ImmuneScore. The x-axis and y-axis are time and survival probability. The blue is High/Low, the red is Low/High, and the green is the Mix group

Full size image

Then the patients were grouped into high-ImmuneScore group (ImmuneScore > 1001.68) and low-ImmuneScore group (ImmuneScore < 1001.68) according to the cutoff value (Fig. 1D). Survival analysis results showed that the survival probability of the low-ImmuneScore group was lower than that of the high-ImmuneScore group (Fig. 1E) (p = 0.059).

Then the patients were reassigned into low-glycolysis-high-ImmuneScore (Low/High) group, high-glycolysis-low-ImmuneScore (High/Low) group and Mix group according to the ImmuneScore and glycolysis. The prognosis of Low/High was better than that of the High/Low group while the prognosis of Mix group was in the middle of the two groups and the prognosis of the three groups differed significantly (Fig. 1F) (p = 0.015).

Screening of glycolysis-immune-related genes

Differential expression analysis was applied to the expression data of OSCC patients to screen the differentially expressed genes (DEGs). There were 2505 DEGs (Fig. 2A) between the low-glycolysis and high-glycolysis groups, 6565 DEGs (Fig. 2B) between low-ImmuneScore group and high-ImmuneScore group and 8503 DEGs (Fig. 2C) between the Low/High and High/Low group. And 337 overlap genes among the three sets of DEGs (Additional file 2: Table S2, Fig. 2D) were found.

Fig. 2
figure 2

Results of differential analysis and enrichment analysis. A–C The volcano plot of DEGs in different groups. The x-axis and y-axis are the multiple of differential expression (log2FC) and − log10(adj. p. val). The blue and orange represent downregulation and upregulation genes. D The figure containing all DEGs. E The significantly enriched GO terms. The x-axis is the number of enriched genes and y-axis is the names of the GO terms. F The significantly enriched KEGG pathways. The x-axis is the number of enriched genes and y-axis is the names of the KEGG pathways

Full size image

Enrichment analysis was performed on the 337 overlap genes. The enriched GO terms and KEGG pathways were listed in Additional file 3: Table S3 (the pathways were obtained basing on KEGG [24,25,26]). The significantly enriched GO terms and KEGG pathways were shown on Fig. 2E and F. The results demonstrated that the 337 overlap genes were significantly enriched in the secondary metabolic process, apical part of cell and carboxylic acid binding GO terms, as well as steroid hormone biosynthesis and HIF-1 signaling KEGG pathways. The genes were significantly enriched in the pathways which were associated with metabolism, inflammation and cancers [27,28,29,30].

Construction and validation of prognostic model

In the TCGA dataset, 337 overlap genes were used as continuous variable to perform univariate Cox regression analysis. Hazard ratio (HR) of each gene was calculated and p < 0.01 was used as threshold to screen the genes which were associated with OS of OSCC (Fig. 3A). Seven genes were selected and then LASSO Cox regression analysis was applied to the 7 genes to screen a set of genes which were significantly associated with the prognosis of OSCC. Six genes (including ALDOC, VEGFA, HRG, PADI3 IGSF11 and MIPOL1) were screened according to the lowest lambda value (Fig. 3B).

Fig. 3
figure 3

Construction of OSCC prognostic model. A The forest plot of seven genes significantly associated with the prognosis of OSCC by univariate Cox regression analysis. HR is short for Hazard ratio and 95%CI is 95% of the confidence interval. B The graph to determine the best lambda by LASSO Cox model. The x-axis is the log (lambda) and y-axis is the partial likelihood deviance whose smallest value is the best lambda value. C The Kaplan–Meier survival curve in TCGA dataset. The x-axis is time and y-axis is the survival probability. The red is the low-risk group and the blue is the high-risk group. p value is calculated by the log-rank test. D The Kaplan–Meier survival curve in meta-GEO dataset. The x-axis is time and y-axis is the survival probability. The red is the low-risk group and the blue is the high-risk group. p value is calculated by the log-rank test. E and F The estimated Risk Score of each patient is ranked from small to big. The vertical dotted line is the median of Risk Score

Full size image

In order to obtain a uniform critical value, we standardized the the 6 gene expression values in the TCGA and three GEO datasets to a value with median 0 and standard deviation 1. Then the standardized expression values were multiplying regression coefficient to construct a Risk Score model as follows: Risk Score = (0.039770315*ALDOC) + (0.164774576*VEGFA) + (0.208429564*HRG) + (0.175124998*PADI3) + (0.065430137*IGSF11) + (0.001152926*MIPOL1). We calculated Risk Score of each patient. And patients in the TCGA and meta-GEO datasets were assigned into high-risk and low-risk groups according to the best cutoff value of Risk Score (− 0.3075). OS of patients in the high-risk group was lower than that in the low-risk group both in the TCGA and meta-GEO datasets (including GSE85446, GSE65858 and GSE41613) from the survival analysis (Fig. 3C–D).

According to the survival time in combination with patients which were ranked by the Risk Score in TCGA and meta-GEO datasets, there were more patients in the high-risk group compared to the low-risk group. In Fig. 3E and F, the green dots and red dots represent alive and dead patients. And the number of dead patients in the high-risk group was higher than that in the low-risk group.

In conclusion, the results suggested that Risk Score model constructed by ALDOC, VEGFA, HRG, PADI3 IGSF11 and MIPOL1 might be a reliable prognostic model.

Risk Score as an independent prognostic hallmark of OSCC

Multivariate Cox regression analysis was then conducted including totally 10 factors, comprising Risk Score, age, gender, Grade, tobacco history, alcohol history, tumor stage, TNM status, to determine the independent prognostic indicators for OSCC patients (Fig. 4A) (removing one non-differentiated sample and 8 samples without stage & T status information). The results showed that the Risk Score, age, Grade, and Node status were significantly correlated with the OS of OSCC (Fig. 4A). Patients in the low-risk group had lower death risk and Risk Score was a reliable prognostic factor (HR = 4.25, 95%CI: 2.676—6.7, p < 0.001).

Fig. 4
figure 4

Risk Score is an independent prognostic biomarker for OSCC. A Forest plot of multiple Cox regression analysis. Compared to the reference, patients with Hazard ratio > 1 were considered higher death risk, and patients with Hazard ratio < 1 were considered lower death risk. B, C Kaplan–Meier survival curves in female and male subgroups, separately. D–E Kaplan–Meier survival curves of ≤ 61-year old and > 61-year old subgroups, separately. F–H Kaplan–Meier survival curves of patients in N0, N1, N2-N3 subgroups, separately

Full size image

To further explore the prognostic values of Risk Score in different pathological factors such as age, gender, and Node status, we regrouped the patients to perform Kaplan–Meier survival analysis. Our data suggested that high risk OSCC patients in various subgroups all had inferior survival compared to the patients in low-risk group, including female patients (Fig. 4B), male patients (Fig. 4C), younger patients (≤ 61-year old) (Fig. 4D), older patients (> 61-year old) (Fig. 4E), and N0, N1, N2-N3 patients (Fig. 4F-H). Above results indicated that the Risk Score was an independent prognostic indicator for stratifying the OSCC patients with different prognosis.

Nomogram model predicts the survival probability of OSCC patients with good performance

Subsequently, the independent prognostic factors, Risk Score, age, Grade, and Node status, were used to construct the nomogram model (Fig. 5A) to predict the survival probability of OSCC patients in 1 year, 3 years and 5 years. In the calibration graph, the adjusted curves were close to the ideal curve (a line through the origin with the slop 1 and 45 degree). These indicated that the estimated survival probability agreed well with the actually living probability (Fig. 5B–D) in 1 year, 3 years and 5 years.

Fig. 5
figure 5

Nomogram model to predict the prognosis of OSCC patients. A Nomogram to the living probability of OSCC patients in 1 year, 3 years and 5 years. B–D The calibration curve of nomogram for OSCC patients in 1 year, 3 years and 5 years, respectively. The x-axis is the survival probability estimated by nomogram and y-axis is the actual survival probability

Full size image

The immunosuppressive cells of OSCC patients were significantly infiltrated in the low-risk group

After organizing the immune infiltrations of 306 OSCC patients (Fig. 6A), we could find that the proportions of the immune infiltration between the high- and low-risk groups differed significantly (Fig. 6B). The infiltration ratios of B cells naive, T cells CD8, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, monocytes, macrophages M0, macrophages M1, macrophages M2, dendritic cells activated, mast cells testing and mast cells activated differed significantly (Fig. 6C). The infiltration ratios of B cells naive, T cells CD8, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, monocytes, macrophages M1, macrophages M2 and mast cells testing were higher in the low-risk group, while the infiltration ratios of macrophages M0, dendritic cells activated and mast cells activated were higher in the high-risk group. The correlation of different immune cells was relatively weak, which might imply the weak interactions among immune cells in OSCC (Fig. 6D).

Fig. 6
figure 6

Immune infiltration between the high- and low-risk groups. A The relative infiltration ratio of immune cells in all OSCC patients. B The violin diagram of immune cells of the difference between the high- and low-risk groups. The x-axis is the 21 immune cells and the y-axis is the infiltration ration of immune cells. We used the wilcoxn method to calculate the p value. (ns: p > 0.05, *: p <  = 0.05, **: p <  = 0.01, ***: p <  = 0.001, ****: p <  = 0.0001). C The box graph of differentially expressed immune cells between the high- and low-risk groups. The x-axis is the high- and low-risk groups. The y-axis is the infiltration ratio of immune cells and p value is calculated by the wilcoxn method. D The correlation matrix of 22 immune cells. The orange is positive correlation and the blue represents negative correlation. The greatness of the correlation is positively related to the darkness of the color. E The chord diagram of Risk Score and 5 immune checkpoints. The width of the connecting lines represent the strength of the correlation. F The violin diagram of the difference of expression of the checkpoints between the high- and low-risk groups. The red is high-risk group and the blue is low-risk group. The y-axis is the expression of the checkpoints, and p value is calculated by the wilcoxn method

Full size image

We analyzed the correlation between the Risk Score and the immune checkpoints (CTLA4, PDL1, LAG3, TIGIT, IDO1 and TDO2), and the results showed that the Risk Score related to 6 immune checkpoints (Fig. 6E). The expressions of LAG3 (p = 9.9*e−5), TIGIT (p = 0.0033) and IDO1 (p = 0.0027) were significantly higher in the low-risk group compared to the high-risk group (Fig. 6F). Moreover, higher PDL1 expression was observed in low-risk OSCC patients, and similar PDL1 expression tendency had also been documented in previous reports [31, 32].

Discussion

OSCC has been widely considered the most frequent malignancy of the head and neck, novel prognostic biomarkers for OSCC are helpful for stratifying patients with different prognosis [33]. In this study, basing on the public OSCC data, glycolysis-immune related genes were employed to build predictive prognostic model for OSCC. Our Risk Score, based on ALDOC, VEGFA, HRG, PADI3, IGSF11 and MIPOL1, was a promising prognostic indicator for OSCC.

Firstly, we downloaded data of OSCC patients from public databases to perform differential expression analysis and functional enrichment analysis. There were 337 overlap DEGs among the groups defined by glycolysis and ImmuneScore, and these genes were enriched in 99 GO terms and 9 KEGG pathways. We noticed that several metabolism related terms were significantly enriched, including secondary metabolic process, apical part of cell and carboxylic acid binding GO terms and the steroid hormone biosynthesis pathway. The carboxylic acid is an important substance in the generation and progress of cancers [34]. Besides, the biosynthesis of steroid hormone, a process requiring multiple enzymes to coordinate [35], is also found to be related to the formation and growth of prostate cancer [36]. Whereas, their roles in OSCC have remained largely unclear. Additionally, the HIF-1 signaling pathway was also significantly enriched. Hypoxia is a common feature of many tumors. Under hypoxia environment, HIF-1 is able to bind to the hypoxia response elements (HREs) of target genes, and the target genes included the genes encoding glycolytic receptors and enzymes [37]. Thus, HIF-1 signaling pathway might imply the indirect correlation between hypoxia and our Risk Score, whose details deserve further investigation.

Then we optimized the DEGs using univariate Cox and LASSO Cox regression analysis and 6 genes were selected, including ALDOC, VEGFA, HRG, PADI3, IGSF11 and MIPOL1. Previous studies have shown that the ALDOC is a member of the aldolase family and has been identified as an independent prognostic hallmark for cancers [38]. Besides, ALDOC has been evidenced to inhibit the migration of OSCC and serve as a prognosis marker [39], which was in line with our data. The VEGFA is responsible for the formation of new blood vessels, which is important for the progression of cancers [40]. VEGFA has also been reported as single gene prognosis marker in OSCC [41], while we firstly included it in the glycolysis-immune related prognostic Risk Score. The HRG expression is connected with prognosis of colorectal cancer [42]. The HRG has been used as prognostic biomarker for many kinds of cancers such as prostate cancer [43]. However, few HRG related studies were found in OSCC. The PADI3 has anticancer effect through the arresting of cell cycle in colon cancer, and it regulates the glycolysis in multiple cancer cell types [44]. The MIPOL1 may induce tumor suppression in nasopharyngeal carcinoma resulting anticancer effect [45]. Although two genes (HRG and PADI3) have been seldom studied in OSCC or HNSCC, their roles in other tumors could be found, indicating that more related exploration should be done in OSCC. The above evidence indicated that these genes, which were associated with glycolysis and immunity, are related to progression of cancers, which suggesting that the glycolysis-immune-related genes might be reliable prognostic biomarkers for OSCC patients.

Our prognostic Risk Score was constructed basing on the 6 core genes to separate patients with different prognosis, and patients in the TCGA and GEO databases were assigned to high- and low-risk groups according to the best cutoff of Risk Score. The Kaplan–Meier survival demonstrated that patients in the high-risk group had lower survival rates than those in the low-risk group, indicating the relatively reliable performance of the Risk Score. The Risk Score models constructed by glycolysis-related genes have been proven prognostic values in pancreatic ductal adenocarcinoma [46]. Prognostic models basing on immune-related genes have also been identified in renal papillary cell carcinoma [47]. The prognostic models constructed by glycolysis-related genes or immune-related genes have been proven reliable in many types of cancers, which are in accord with our results. Additionally, between high and low risk OSCC patients, 12 types of immune cells were significantly differentially infiltrated, indicating the different immune microenvironment between the patients with different prognosis. Herein, relatively higher abundance of multiple immune cells were observed in low risk patients, including Treg and Macrophages M2 cells. Tregs have been reported to inhibit the antigen-presenting cells and thereby promote the proliferation of tumor cells, meanwhile Tregs are also related to undesirable prognosis [48]. However, more recently, it has been indicated that the infiltration, activation, and survival of Tregs involve in complicated multiple processes in TME, which differ among distinct tumor types [48]. Moreover, complicated factors, including the composition and activity of infiltrated immune cells in tumor immune microenvironment and the cell surface expression of immune checkpoints, together determine the immune response states in tumor microenvironment [49, 50] Several key immune checkpoints also showed distinct expression between high and low risk OSCC patients. To the best of our knowledge, this is the first study assessing the prognostic value of glycolysis-immune-related genes for OSCC patients, known to be implicated in stratifying patients with different prognosis.

Nevertheless, some limitations of our study were needed to be further improved in the near future. Although we have validated our prognostic Risk Score in TCGA and meta-GEO datasets, more validation in the expanded sample size might further elevate the confidence of the Risk Score. Additionally, more underlying functional details about the 6 core genes of Risk Score in OSCC should be investigated in our future work.

Conclusions

In this work, we have firstly revealed a reliable glycolysis-immune related prognostic Risk Score for OSCC patients, basing on multiple OSCC public datasets. The prognostic Risk Score was based on 6 core genes, comprising ALDOC, VEGFA, HRG, PADI3, IGSF11 and MIPOL1. Our prognostic Risk Score model might be helpful in separating OSCC patients with different prognosis, shedding new light on the future management choosing of OSCC patients.

Availability of data and materials

The datasets analysed during the current study are available in The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/) and the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/, accession number: GSE85446, GSE65858 and GSE41613) databases.

Abbreviations

ATP:

Adenosine triphosphate

DEGs:

Differentially expressed genes

GEO:

Gene expression omnibus

GO:

Gene ontology

HR:

Hazard ratio

HPV:

Human papillomavirus

LDHA:

Lactate dehydrogenase A

KEGG:

Kyoto encyclopedia of genes and genomes

OSCC:

ORAL squamous cell carcinoma

OS:

Overall survival

SSE:

Sum of the square errors

TCGA:

The cancer genome atlas

References

  1. Chamoli A, Gosavi AS, Shirwadkar UP, Wangdale KV, Behera SK, Kurrey NK, et al. Overview of oral cavity squamous cell carcinoma: risk factors, mechanisms, and diagnostics. Oral Oncol. 2021;121:105451.

    Article  PubMed  Google Scholar 

  2. McCord C, Kiss A, Magalhaes MA, Leong IT, Jorden T, Bradley G. Oral squamous cell carcinoma associated with precursor lesions. Cancer Prev Res. 2021;14(9):873–84 (Phila).

    Article  Google Scholar 

  3. Ferreira AK, Carvalho SH, Granville-Garcia AF, Sarmento DJ, Agripino GG, Abreu MH, et al. Survival and prognostic factors in patients with oral squamous cell carcinoma. Med Oral Patol Oral Cir Bucal. 2021;26(3):e387–92.

    Article  PubMed  Google Scholar 

  4. Panarese I, Aquino G, Ronchi A, Longo F, Montella M, Cozzolino I, et al. Oral and Oropharyngeal squamous cell carcinoma: prognostic and predictive parameters in the etiopathogenetic route. Expert Rev Anticancer Ther. 2019;19(2):105–19.

    Article  PubMed  Google Scholar 

  5. Yang Z, Yan G, Zheng L, Gu W, Liu F, Chen W, et al. YKT6, as a potential predictor of prognosis and immunotherapy response for oral squamous cell carcinoma, is related to cell invasion, metastasis, and CD8+ T cell infiltration. Oncoimmunology. 2021;10(1):1938890.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Russo D, Merolla F, Mascolo M, Ilardi G, Romano S, Varricchio S, et al. FKBP51 immunohistochemical expression: a new prognostic biomarker for OSCC? Int J Mol Sci. 2017;18(2):443.

    Article  PubMed Central  Google Scholar 

  7. Dienel GA. Brain glucose metabolism: integration of energetics with function. Physiol Rev. 2019;99(1):949–1045.

    Article  PubMed  Google Scholar 

  8. Chaudhry R, Varacallo M. Biochemistry, Glycolysis. StatPearls. Treasure Island (FL). 2022.

  9. Melkonian EA, Schury MP. Biochemistry, Anaerobic Glycolysis. StatPearls. Treasure Island (FL). 2022.

  10. Vaupel P, Schmidberger H, Mayer A. The Warburg effect: essential part of metabolic reprogramming and central contributor to cancer progression. Int J Radiat Biol. 2019;95(7):912–9.

    Article  PubMed  Google Scholar 

  11. Eckert AW, Lautner MH, Schutze A, Taubert H, Schubert J, Bilkenroth U. Coexpression of hypoxia-inducible factor-1alpha and glucose transporter-1 is associated with poor prognosis in oral squamous cell carcinoma patients. Histopathology. 2011;58(7):1136–47.

    Article  PubMed  Google Scholar 

  12. Gong X, Tang H, Yang K. PER1 suppresses glycolysis and cell proliferation in oral squamous cell carcinoma via the PER1/RACK1/PI3K signaling complex. Cell Death Dis. 2021;12(3):276.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Wang Y, Li Y, Jiang L, Ren X, Cheng B, Xia J. Prognostic value of glycolysis markers in head and neck squamous cell carcinoma: a meta-analysis. Aging. 2021;13(5):7284–99 (Albany NY).

    Article  PubMed  PubMed Central  Google Scholar 

  14. Li Z, Zheng J, Feng Y, Li Y, Liang Y, Liu Y, et al. Integrated analysis identifies a novel lncRNA prognostic signature associated with aerobic glycolysis and hub pathways in breast cancer. Cancer Med. 2021;10(21):7877–92.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Chang CH, Curtis JD, Maggi LB Jr, Faubert B, Villarino AV, O’Sullivan D, et al. Posttranscriptional control of T cell effector function by aerobic glycolysis. Cell. 2013;153(6):1239–51.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Xu K, Yin N, Peng M, Stamatiades EG, Shyu A, Li P, et al. Glycolysis fuels phosphoinositide 3-kinase signaling to bolster T cell immunity. Science. 2021;371(6527):405–10.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Brand A, Singer K, Koehl GE, Kolitzus M, Schoenhammer G, Thiel A, et al. LDHA-associated lactic acid production blunts tumor immunosurveillance by T and NK cells. Cell Metab. 2016;24(5):657–71.

    Article  PubMed  Google Scholar 

  18. Keulers AR, Kiesow L, Mahnken AH. Port implantation in patients with severe thrombocytopenia is safe with interventional radiology. Cardiovasc Intervent Radiol. 2018;41(1):80–6.

    Article  PubMed  Google Scholar 

  19. Peng L, Chen Z, Chen Y, Wang X, Tang N. MIR155HG is a prognostic biomarker and associated with immune infiltration and immune checkpoint molecules expression in multiple cancers. Cancer Med. 2019;8(17):7161–73.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51.

    Article  PubMed  Google Scholar 

  27. Balamurugan K. HIF-1 at the crossroads of hypoxia, inflammation, and cancer. Int J Cancer. 2016;138(5):1058–66.

    Article  PubMed  Google Scholar 

  28. Manikandan P, Nagini S. Cytochrome P450 structure, function and clinical significance: a review. Curr Drug Targets. 2018;19(1):38–54.

    Article  PubMed  Google Scholar 

  29. Croom E. Metabolism of xenobiotics of human environments. Prog Mol Biol Transl Sci. 2012;112:31–88.

    Article  PubMed  Google Scholar 

  30. Ediriweera MK, Tennekoon KH, Samarakoon SR. Role of the PI3K/AKT/mTOR signaling pathway in ovarian cancer: biological and therapeutic significance. Semin Cancer Biol. 2019;59:147–60.

    Article  PubMed  Google Scholar 

  31. Xu F, Zhan X, Zheng X, Xu H, Li Y, Huang X, et al. A signature of immune-related gene pairs predicts oncologic outcomes and response to immunotherapy in lung adenocarcinoma. Genomics. 2020;112(6):4675–83.

    Article  PubMed  Google Scholar 

  32. Jiang P, Li Y, Xu Z, He S. A signature of 17 immune-related gene pairs predicts prognosis and immune status in HNSCC patients. Transl Oncol. 2021;14(1):100924.

    Article  PubMed  Google Scholar 

  33. De Paz D, Kao HK, Huang Y, Chang KP. Prognostic stratification of patients with advanced oral cavity squamous cell carcinoma. Curr Oncol Rep. 2017;19(10):65.

    Article  PubMed  Google Scholar 

  34. Bian X, Qian Y, Tan B, Li K, Hong X, Wong CC, et al. In-depth mapping carboxylic acid metabolome reveals the potential biomarkers in colorectal cancer through characteristic fragment ions and metabolic flux. Anal Chim Acta. 2020;1128:62–71.

    Article  PubMed  Google Scholar 

  35. Bacila IA, Elder C, Krone N. Update on adrenal steroid hormone biosynthesis and clinical implications. Arch Dis Child. 2019;104(12):1223–8.

    Article  PubMed  Google Scholar 

  36. Neuwirt H, Bouchal J, Kharaishvili G, Ploner C, Johrer K, Pitterl F, et al. Cancer-associated fibroblasts promote prostate tumor growth and progression through upregulation of cholesterol and steroid biosynthesis. Cell Commun Signal. 2020;18(1):11.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Zheng F, Chen J, Zhang X, Wang Z, Chen J, Lin X, et al. The HIF-1alpha antisense long non-coding RNA drives a positive feedback loop of HIF-1alpha mediated transactivation and glycolysis. Nat Commun. 2021;12(1):1341.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Chang YC, Yang YC, Tien CP, Yang CJ, Hsiao M. Roles of aldolase family genes in human cancers and diseases. Trends Endocrinol Metab. 2018;29(8):549–59.

    Article  PubMed  Google Scholar 

  39. Li YJ, Huang TH, Hsiao M, Lin BR, Cheng SJ, Yang CN, et al. Suppression of fructose-bisphosphate aldolase C expression as a predictor of advanced oral squamous cell carcinoma. Head Neck. 2016;38(Suppl 1):E1075–85.

    Article  PubMed  Google Scholar 

  40. Cuzziol CI, Castanhole-Nunes MMU, Pavarino EC, Goloni-Bertollo EM. MicroRNAs as regulators of VEGFA and NFE2L2 in cancer. Gene. 2020;759:144994.

    Article  PubMed  Google Scholar 

  41. Peterle GT, Maia LL, Trivilin LO, de Oliveira MM, Dos Santos JG, Mendes SO, et al. PAI-1, CAIX, and VEGFA expressions as prognosis markers in oral squamous cell carcinoma. J Oral Pathol Med. 2018;47(6):566–74.

    Article  PubMed  Google Scholar 

  42. Lee JH, Jung S, Park WS, Choe EK, Kim E, Shin R, et al. Prognostic nomogram of hypoxia-related genes predicting overall survival of colorectal cancer-analysis of TCGA database. Sci Rep. 2019;9(1):1803.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Yang L, Roberts D, Takhar M, Erho N, Bibby BAS, Thiruthaneeswaran N, et al. Development and validation of a 28-gene hypoxia-related prognostic signature for localized prostate cancer. EBioMedicine. 2018;31:182–9.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Coassolo S, Davidson G, Negroni L, Gambi G, Daujat S, Romier C, et al. Citrullination of pyruvate kinase M2 by PADI1 and PADI3 regulates glycolysis and cancer cell proliferation. Nat Commun. 2021;12(1):1718.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Leong MML, Cheung AKL, Kwok TCT, Lung ML. Functional characterization of a candidate tumor suppressor gene, mirror image polydactyly 1, in nasopharyngeal carcinoma. Int J Cancer. 2020;146(10):2891–900.

    Article  PubMed  Google Scholar 

  46. Song W, He X, Gong P, Yang Y, Huang S, Zeng Y, et al. Glycolysis-related gene expression profiling screen for prognostic risk signature of pancreatic ductal adenocarcinoma. Front Genet. 2021;12:639246.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Wang Z, Song Q, Yang Z, Chen J, Shang J, Ju W. Construction of immune-related risk signature for renal papillary cell carcinoma. Cancer Med. 2019;8(1):289–304.

    Article  PubMed  Google Scholar 

  48. Nishikawa H, Koyama S. Mechanisms of regulatory T cell infiltration in tumors: implications for innovative immune precision therapies. J Immunother Cancer. 2021;9(7):e002591.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Zhang Y, Liu Q, Liao Q. Long noncoding RNA: a dazzling dancer in tumor immune microenvironment. J Exp Clin Cancer Res. 2020;39(1):231.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Mao X, Xu J, Wang W, Liang C, Hua J, Liu J, et al. Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Cancer. 2021;20(1):131.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Author information

Authors and Affiliations

  1. Department of Stomatology, Tianjin First Central Hospital, Nankai District, No.24 Fukang Road, Tianjin, 300192, People’s Republic of China

    Yi Liu, Tong Wang & Ronghua Li

Authors
  1. Yi Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ronghua Li
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YL contributed to the study conception and design. Data collection and analysis were performed by all authors. The first draft of the manuscript was written by YL and TW. RHL commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Yi Liu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

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

Supplementary Information

Additional file 1

. Table S1. The glycolysis-related genes from GSEA database.

Additional file 2

. Table S2. The overlap genes of the DEGs screened by levels of glycolysis and Immune Score.

Additional file 3

. Table S3. The enriched GO terms and KEGG pathways.

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, Y., Wang, T. & Li, R. A prognostic Risk Score model for oral squamous cell carcinoma constructed by 6 glycolysis-immune-related genes. BMC Oral Health 22, 324 (2022). https://doi.org/10.1186/s12903-022-02358-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12903-022-02358-0

Keywords

BMC Oral Health

ISSN: 1472-6831

Contact us