- Research
- Open access
- Published:
Construction and validation of a metabolic-related genes prognostic model for oral squamous cell carcinoma based on bioinformatics
BMC Medical Genomics volume 15, Article number: 269 (2022)
Abstract
Background
Oral squamous cell carcinoma (OSCC) accounts for a frequently-occurring head and neck cancer, which is characterized by high rates of morbidity and mortality. Metabolism-related genes (MRGs) show close association with OSCC development, metastasis and progression, so we constructed an MRGs-based OSCC prognosis model for evaluating OSCC prognostic outcome.
Methods
This work obtained gene expression profile as well as the relevant clinical information from the The Cancer Genome Atlas (TCGA) database, determined the MRGs related to OSCC by difference analysis, screened the prognosis-related MRGs by performing univariate Cox analysis, and used such identified MRGs for constructing the OSCC prognosis prediction model through Lasso-Cox regression. Besides, we validated the model with the GSE41613 dataset based on Gene Expression Omnibus (GEO) database.
Results
The present work screened 317 differentially expressed MRGs from the database, identified 12 OSCC prognostic MRGs through univariate Cox regression, and then established a clinical prognostic model composed of 11 MRGs by Lasso-Cox analysis. Based on the optimal risk score threshold, cases were classified as low- or high-risk group. As suggested by Kaplan–Meier (KM) analysis, survival rate was obviously different between the two groups in the TCGA training set (P < 0.001). According to subsequent univariate and multivariate Cox regression, risk score served as the factor to predict prognosis relative to additional clinical features (P < 0.001). Besides, area under ROC curve (AUC) values for patient survival at 1, 3 and 5 years were determined as 0.63, 0.70, and 0.76, separately, indicating that the prognostic model has good predictive accuracy. Then, we validated this clinical prognostic model using GSE41613. To enhance our model prediction accuracy, age, gender, risk score together with TNM stage were incorporated in a nomogram. As indicated by results of ROC curve and calibration curve analyses, the as-constructed nomogram had enhanced prediction accuracy compared with clinicopathological features alone, besides, combining clinicopathological characteristics with risk score contributed to predicting patient prognosis and guiding clinical decision-making.
Conclusion
In this study, 11 MRGs prognostic models based on TCGA database showed superior predictive performance and had a certain clinical application prospect in guiding individualized.
Background
Oral cancer (OC) accounts for the head and neck cancer that mainly originates from the cheek, tongue, lip and palate. There are more than 300,000 new oral cancer cases and more than 145,000 related deaths annually [1]. Oral squamous cell carcinoma (OSCC) occupies ninety percent of OC patients. The traditional treatment of OSCC is mainly surgical resection of the tumor, preoperative or postoperative chemotherapy, radiotherapy and adjuvant therapy [2]. As researchers continue to study this disease, the treatment options for these patients continue to improve. Image-based and adaptive radiotherapy, transoral robotic resection and immunotherapy are gradually being utilized to treat OSCC [3,4,5]. Unfortunately, most patients are diagnosed at a late stage, the local recurrence rate is high, and metastasis often occurs, making the 5-year survival rate much lower than that of other malignant tumors. In fact, during the previous 30 years, the survival rate of oral squamous cell carcinoma patients has been consistently less than 50% [6, 7]. More accurate prediction of the prognosis of OSCC patients will allow doctors to better choose appropriate treatment strategies and improve the survival rate of patients [8]. Currently, the tumor, node, metastasis (TNM) classification system has been mainly used to predict tumor prognosis and provide clinical guidance in choosing appropriate treatment methods. However, the OSCC-related TNM system mainly focuses on the anatomical extent of the disease and ignores factors related to tumor prognosis, such as age, gender, and the presence of other diseases. Thus, patients who have the same TNM stage can have very different survival outcomes; in other words, the use of the TNM staging alone for the prediction of patient survival is insufficient [9]. The eighth American Joint Commission on Cancer (AJCC) staging system emphasizes the need for a "personalized" treatment approach for cancer patients [10]. In addition, inaccurate prognostic information can affect treatment decisions and subsequent outcomes. For example, high-risk cases are possibly associated with cancer migration and metastasis due to insufficient or delayed treatment, and low-risk patients may experience a loss of bone marrow function and organ function because of excessive treatment; both situations have substantial effect on patient treatment or recovery [11]. Consequently, it is urgently needed to construct the creditable prognostic approaches that can assist clinicians in selecting the suitable individualized therapeutic strategies, thus improving OSCC prognosis.
Changes in metabolic processes are closely related to tumor growth [12,13,14,15]. During the process of rapid growth and proliferation of tumor cells, the metabolic pathways of the body are in constant flux. In other words, during the process of tumorigenesis, progression and metastasis, tumor cells reprogram catabolic and anabolic pathways to satisfy the energy metabolism and biosynthesis of cells by enhancing macromolecular biosynthesis, regulating the redox balance and allowing for the rapid production of ATP [16, 17]. In the 1920s, Otto Warburg discovered that under both aerobic and anaerobic conditions, the ATP produced by tumor cells tends to undergo glycolysis rather than oxidative phosphorylation, which offers more energy for the rapid proliferation of tumor cells; this has been termed the “Warburg effect” [18]. In addition, the metabolic processes involved in tumor metabolic reprogramming include the pentose phosphate pathway, lipid metabolism, and nucleic acid and amino acid metabolism [19]. During tumorigenesis, progression and metastasis, changes in MRGs determine the changes in metabolic pathways [20]. For example, some studies have shown that PER1 and PER2 are significantly downregulated in OSCC and that their overexpression can inhibit glycolysis and tumor cell proliferation, thus inhibiting tumor growth [21, 22]. In contrast to PER1 and PER2, the expression of SHMT2 in OSCC was found to be significantly upregulated, which predicts the dismal OSCC survival [23]. In addition, the upregulation of PFKP in OSCC is related to the pathological differentiation of tumors as well as lymph node metastasis [24]. According to the obtained results, metabolism-related genes (MRGs) are a prognostic marker and therapeutic target for OSCC.
As bioinformatics analysis is increasingly used for diagnosing and predicting cancer prognostic outcome, some researchers have connected metabolome with genome [25]. To date, many studies have used MRGs to construct clinical prognostic models of malignant tumors like gastric cancer, cervical cancer, liver cancer and renal clear cell carcinoma, and these studies have achieved good predictive results [26,27,28]. However, the development of clinical prognostic models for OSCC based on MRGs has been barely explored and are not comprehensive [29]. In this study, clinical data and gene levels for OSCC cases were acquired in The Cancer Genome Atlas (TCGA). This TCGA-OSCC cohort was then used as the training set to screen the differentially expressed MRGs that were closely related to OSCC by difference analysis. The MRGs found to be significantly related to prognosis were screened by univariate Cox analysis, and the clinical prognosis model of OSCC based on 11 MRGs was constructed by Lasso-Cox analysis. Finally, risk score was calculated and identified as the factor that independently predicted OSCC prognosis upon univariate as well as multivariate Cox regression. In addition, we further validated the dataset obtained using the Gene Expression Omnibus (GEO) database as the validation set. This study describes a new approach to predict overall survival (OS) of OSCC cases, which could help direct the individualized treatment of OSCC patients.
Materials and methods
Data downloading and processing
This work obtained clinical data together with the FPKM-normalized gene profiles for 340 cancers along with 32 non-carcinoma samples in TCGA database OSCC cohort (https://portal.gdc.cancer.gov/). Apart from that, this work also obtained the GSE41613 dataset, which includes 97 OSCC samples, in GEO database (https://www.ncbi.nlm.nih.gov/geo/). According to annotation patterns obtained by relevant platform, this work transformed the gene matrix file of probe identification as gene symbols. In addition, samples collected in cases with a < 90-day follow-up period from TCGA-OSCC cohort and GSE41613 were excluded, and the survival time was changed from years/months to days. The TCGA-OSCC cohort was enrolled into training cohort, while the GSE41613 dataset into external validation cohort. This work acquired a total of 851 MRGs (c2.cp.kegg.v7.4.symbols.gmt) in GSEA (http://www.gsea-msigdb.org/gsea/index.jsp).
Screening differentially expressed genes (DEGs)
After genes from training cohort were intersected with those from metabolic gene cohort, the MRGs in the TCGA gene expression profile were extracted. This work employed R software (Version 4.1.2) limma R package for screening metabolic DEGs within cancer compared with non-carcinoma samples upon |log2FC|> 0.5 and FDR < 0.05 thresholds. Then, we used the glmnet R package to integrate survival status, survival time, as well as DEGs expression profiles.
Construction of the MRGs-based prognostic signature
By univariate Cox regression, this work calculated hazard ratios (HRs) as well as relevant 95% confidence intervals (CIs) and obtained 12 genes related to overall survival. Next, we selected the most appropriate lambda value to obtain the optimal model through lasso analysis, and screened the above 12 genes into 11 MRGs. Meanwhile, the present study calculated risk score by the gene correlation coefficient obtained by Lasso-Cox analysis. Survminer R and survival R packages were utilized to conduct Kaplan–Meier (KM) analysis, and the cutoff value of continuous variables in survival data was measured through survminer R package surv_cutpoint function. In addition, pROC package roc function was also adopted for analyzing 1, 3 and 5 years receiver operating characteristic (ROC) curves, separately, whereas area under the curve (AUC) values and CIs were adopted for calculating pROC package ci function for obtaining the final AUC results and determining risk score’s sensitivity and accuracy in the prediction of OSCC prognosis at 1, 3 and 5 years. Then, this work carried out univariate as well as multivariate Cox analysis for comparing the association of risk score with additional clinicopathological factors in prognosis prediction for training set (including age, gender, stage, grade, T stage and N stage). Finally, association of risk score based on MRGs with clinical factors was analyzed for better confirming that our prognosis model was accurate. M stage was excluded because of the lack of data.
Taking GSE41613 as the verification set, this work determined risk score by the calculation formula used for training cohort. All cases were classified as low- and high-risk groups based on optimal threshold. ROC survival analysis was carried out on the verification set.
Construction and verification of nomogram
According to age, gender, risk score and TNM stage, R software rms package was adopted to integrate the data of seven characteristics, such as survival status or survival time. This work then predicted 1, 3 and 5 year survival for OSCC cases by adding the total scores of the points of each factor into that as-constructed nomogram. The high score indicates the low survival probability. Thereafter, we evaluated the nomogram performance in prognosis prediction by ROC curves and calibration curves.
Enrichment analysis
With the purpose of understanding the gene function in the model, R software was employed to conduct enrichment analysis by adopting Gene Ontology (GO) along with Kyoto Gene and Genomic Encyclopedia (KEGG) databases to obtain the gene set enrichment results. Then, we used GSEA (Version 4.1.0) software for the analysis. P < 0.05 was defined as a significantly enriched pathway to explore the differences in the low-risk and high-risk groups in the training set.
Immune cell infiltration
Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the relative ratio of infiltrating immune cells.
Single-cell sequencing and cell type identification
ScRNA-Seq data from the GSE172577 was obtained from the GEO database. Seurat package is used to process the expression matrix of. Cells that meet the following conditions (> 15,000 UMI/cell, < 6000 genes/cell, > 6000 genes/cell and > 20% mitochondrial genes) were considered low-quality cells and were removed. DoubletFinder package was used to identify and remove potential doublets. We use principal component analysis (PCA) to reduce the dimensionality of the scRNA-Seq dataset. Top 30 principal components (PCs) were used for UMAP. The main cell clusters were identified with the FindClusters function offered by Seurat.
Statistical analysis
R software (Version 4.1.2) was employed for performing statistical analysis. KM curve analysis and log-rank test were performed for analyzing different survival rates between both groups. p-value was determined according to two-sided statistical tests. Except for the statistical standards specified, all statistically significant results needed to satisfy P < 0.05.
Results
Identification of differentially expressed MRGs
In this study, the training set included 319 patients diagnosed with OSCC, and the GSE41613 validation set included 94 patients diagnosed with OSCC. This work enrolled altogether 413 OSCC cases in both datasets. In the Limma software package, 317 MRGs heatmaps (Fig. 1a) and volcano maps (Fig. 1b) were screened by difference analysis, including 175 up-regulated genes and 142 down-regulated genes.
Construction and verification of the MRGs-based prognostic signature
The expression of each index within training cohort was examined through univariate Cox regression, as shown in Fig. 2a: Twelve differentially expressed MRGs related to prognosis were detected as potential prognostic molecular markers. Red is used to depict the genes that are positively correlated with a poor prognosis, and the results indicate that all 12 prognosis-related genes screened are risk genes (HR > 1). Eleven prognostic gene models were established by Lasso-Cox analysis, which included SHMT2, HPRT1, POLD2, HADHB, POLE3, ADK, GOT1, ATIC, MGST1, ADA and GNPDA1 (Fig. 2b and c). Table 1 lists the correlation coefficients of diverse genes. Risk scores were determined by using the correlation coefficient of the gene, as follows: (0.004 * expression of SHMT2 + 0.003 * expression level of HPRT1 + 0.002 * expression level of POLD2 + 0.016 * expression level of HADHB + 0.002 * expression level of POLE3 + 0.002 * expression level of ADK + 0.001 * expression level of GOT1 + 0.006 * expression level of ATIC + 0.006 * expression level of MGST1 + 0.014 * expression level of ADA + 0.039 * expression level of GNPDA1). For training cohort, its risk scores were 0.566–2.989, then cases were classified as low-risk (n = 160) or high-risk (n = 159) group based on the optimal threshold. KM curve analysis in Fig. 3a was conducted for analyzing survival of low- and high-risk patients. As a result, high-risk patients had markedly decreased survival in comparison with low-risk patients (P < 0.0001). With regard to risk score, expression pattern and survival status distributions (Fig. 3b), from left to right, the risk scores of the patients were sorted from the lowest to the highest, and the dots represent the OSCC patients. Meanwhile, in heatmap, the green-to-red color stands for low-to-high expression level. Based on the above results, the 11 MRGs showed increased levels with risk score, while OSCC survival declined. As revealed by ROC curves at 1, 3 and 5 years (Fig. 3c), the AUC values of training cohorts were determined to be 0.63, 0.70, and 0.76, separately, which indicated that our as-constructed 11 MRGs-based model served as the favorable prognostic tool with good predictive accuracy and sensitivity. We also constructed a ROC curve (Fig. 3d) involving clinical risk factors like age and gender as well as the risk score. The results showed that compared with age (AUC = 0.575), gender (AUC = 0.563), stage (AUC = 0.556), grade (AUC = 0.557), T stage (AUC = 0.559) and N stage (AUC = 0.567), the risk score had a larger AUC (AUC = 0.670), suggesting that the risk score based on MRGs could have a certain role in predicting prognosis.
Samples that had insufficient clinical data and the M stage, and those also had a large amount of missing data (n = 168) were excluded, the clinical information, pathological features and risk scores were analyzed through univariate as well as multivariate Cox regression. As revealed by univariate Cox regression (Fig. 4a), risk score showed significant relation to patient prognosis (HR = 5.125, 95% CI 3.162–8.392, P < 0.001). Besides, clinicopathological characteristics like age, grade and stage were adjusted, as a result, risk score remained the factor independently predicting OSCC upon multivariate Cox regression (Fig. 4b) (P < 0.001). Clinical correlation analysis further compared the relation of risk score with clinicopathological features like age, gender, stage, and grade. Figure 4c shows that the risk scores based on gender, stage and grade in the training set were statistically significant; the risk score of stage III-IV was obviously higher than that of stage I-II (P < 0.05), that of G3-G4 markedly increased compared with G1-G2 (P < 0.05), while that of male patients remarkably increased relative to female patients (P < 0.05). Therefore, the MRGs-based risk score was tightly associated with clinical features, especially staging, grading and gender.
For verifying that our as-constructed model was accurate in predicting OSCC prognosis, we used GSE41613 to verify prognosis model. Cases from verification cohorts were classified as low- (n = 37) or high-risk (n = 57) group based on optimal threshold, by adopting identical formula to that in training cohort to calculate risk score. As suggested by the survival curve (Fig. 5a), high-risk cases had decreased survival compared with low-risk cases (P < 0.001). Figure 5b presents risk score, survival time and survival status distributions between both groups. As a result, the same trend as the training set was observed; that is, an increased risk score predicted the lower survival rate, and the more frequently the patients expressed risk genes. For ROC curves at 1, 3 and 5 years (Fig. 5c), their AUC values were determined to be 0.78, 0.70, and 0.68, separately, while that for ROC curve of risk score (Fig. 5d) was 0.777, indicating that our model better predicted OSCC prognosis than that based on clinical factors alone. Upon univariate as well as multivariate Cox regression (Fig. 6a and b), risk score was significantly correlated with prognosis. The above results are basically consistent with the verification set, which more full (P < 0.001) proves that our model was accurate in prediction.
Construction and verification of nomogram
Based on the above predictors, we constructed a nomogram (Fig. 7). The calibration curves at 1, 3 and 5 years (Fig. 8a–c) approached 45°, which proves that the nomogram is accurate in predicting OSCC survival at 1, 3 and 5 years ROC curve analysis (Fig. 8d–f) was used to test the accuracy of the nomogram score, risk score, and clinical and pathological features at 1, 3 and 5 years. The results showed that relative to additional factors, like age, risk score, T, N stage and grade, the nomogram had higher accuracy in predicting prognosis at 1 year (AUC = 0.710), 3 years (AUC = 0.741) and 5 years(AUC = 0.753). According to the above results, our constructed nomogram well predicted OSCC survival.
Enrichment analysis
The function of 11 MRGs in OSCC biology was explored by GO as well as KEGG analysis. According to GO analysis results (Fig. 9a), those MRGs were mainly linked with purine ribonucleoside monophosphate metabolic process, purine nucleoside monophosphate metabolic process, nucleoside metabolic process, glycosyl compound metabolic process and so on. Upon KEGG analysis (Fig. 9b), MRGs were mostly associated with metabolic pathways, purine metabolism, drug metabolism-other enzymes, carbon metabolism, phenylalanine metabolism and so on. Next, we used GSEA software (Fig. 9c) to select five pathways markedly associated with high-risk patients, such as "β-Alanine metabolism", "cysteine and methionine metabolism", "purine metabolism", and "pyrimidine metabolism". The low-risk group showed α-linolenic acid metabolic enrichment. The above results illustrate the significantly different biological processes in low- versus high-risk patients.
Association between risk signature and immune cell
As shown in Fig. 10, dendritic cells (DC), immature dendritic Cells (iDC), mast cells, and T helper cell 17 (Th17) were more enriched in low-risk group (P < 0.05). T helper cell 2 (Th 2) were more enriched in high-risk group (P < 0.05). DC are the most potent antigen precursor cells in the immune system. DC-mediated cross-initiation of tumor-specific T cells plays a crucial role in initiating and maintaining anti-tumor immunity. Their presence in tumors tends to induce t-cell responses that slow cancer progression and is associated with improved patient survival [30]. Mast cells also play a multifaceted role in the tumor microenvironment by regulating multiple tumor biological events such as cell proliferation and survival, angiogenesis, invasion and metastasis [31]. These findings suggest that immune function is more active in low-risk group. Low-risk patients tend to have a better prognosis.
A single-cell transcriptome atlas in OSCC
31,719 single cells were clustered into seven major cell types through marker genes: epithelial cells (marked with EpCAM, KRT18 and KRT8); T cells (marked with CD2, CD3G, CD3D and CD3E); myeloid cell (marked with LYZ, MS4A6, PECAM1, ENG); fibroblasts (marked with COL6A1, DCN, COL3A1, COL1A1 and COL1A2); endotheliocyte cell (marked with PECAM1, ENG and VWF); B cell (marked with MS4A1, CD79A and CD79B); mast cell (marked with CPA3, KIT and TPSAB1) (Fig. 11a and c). Then, we explored the expression of prognosis-related MRGs in various cells of OSCC (Fig. 11b). We found that five genes are basically expressed in epithelial cells. In other words, as epithelial-derived tumors, the abnormal metabolic patterns in OSCC tumor tissues may be limited to tumor cells. The other six MRGs that make up the prognostic model are not expressed in tumor microenvironment and tumor cells. The reason may be that single cell RNA-seq data usually contain many missing values caused by the failure of original RNA amplification. The proportion of each cell lineage varies greatly among different samples (Fig. 11d).
Disscussion
Oral squamous cell carcinoma originates from oral keratinocytes. It has an insidious onset, making diagnosis difficult, and exhibits rapid development. Due to these characteristics, it is often not detected at an early stage; and the diagnosis of OSCC is usually made at the advanced stage, and the diagnosed patients usually develop distant metastasis upon diagnosis [32, 33]. Thus, accurate prognostic prediction is highly valuable, as it can aid in selecting the most suitable treatment for patients with OSCC, thereby improving their survival [34]. Therefore, it is significant to identify prognosis-related molecular markers comprehensively reflecting tumor biological features, so as to prevent and treat OSCC. With the development of bioinformatics, several large cancer databases, such as TCGA and GEO, have provided researchers access to large-scale gene expression data and corresponding clinical information [35]. Therefore, scholars can explore and develop new, more accurate prognostic models from the perspective of tumor cell biological behavior to optimize the treatment strategy of OSCC. An OSCC prognostic model constructed by using ferroptosis- and autophagy-related genes has been established to investigate the ability in predicting OSCC prognosis and the possibility of targeted treatment [36,37,38]. Many articles suggest the involvement of metabolism in OSCC genesis and progression, but little research regarding MRGs’effect on prognosis prediction of OSCC is available [39,40,41].
To analyze MRGs associated with OSCC survival, this work selected 12 prognostic MRGs from TCGA-OSCC cohort through univariate Cox and differential analyses. A clinical prognostic model of OSCC based on 11 MRGs was established by Lasso-Cox analysis, including 11 genes (SHMT2, HPRT1, POLD2, HADHB, POLE3, ADK, GOT1, ATIC, MGST1, ADA and GNPDA1).
Among these 11 genes, SHMT2, encoded by SHMT2, is an important enzyme related to carbon metabolism in OSCC, as it catalyzes the conversion of serine to glycine [42]. SHMT2 level is discovered to be up-regulated in OSCC tissues, which has been linked to a poor prognosis; the increased expression predicts the worse the pathological state of the tumor [43]. In addition, the silencing of SHMT2 in OSCC cells can affect cell cycle regulatory factors and induce cell G1 phase arrest, resulting in decreased cancer cell growth as well as inhibited cancer proliferation in vivo [44]. SHMT2 is also highly expressed in thyroid cancer, bladder cancer and intrahepatic cholangiocarcinoma, and is closely related to the poor prognosis of these three cancers [45, 46]. HPRT1 participates in regulating cell cycle mainly by regulating the production of purine and inosine in the remedial synthesis pathway [47]. The overexpression of HPRT1 predicts the dismal survival of OSCC and enhances the resistance to cisplatin by promoting the MMP1/PI3K/AKT signaling pathway [48]. In addition, HPRT1 overexpression dramatically decrease immunocyte activities, thus promoting the formation of an immunosuppressive tumor microenvironment [49]. ADA and ADK are mainly involved in adenosine metabolism. In the process of adenosine metabolism, ADA deaminates adenosine to yield inosine, and ADK phosphorylates adenosine to yield adenosine monophosphate [50, 51]. As indicated in this work, the levels of saliva ADA and serum ADA in OSCC cases markedly increased compared with healthy subjects, and the level of serum ADA increased with increasing histopathological grade [52, 53]. Some studies have confirmed that ADA is down-regulated in lymphocytes of advanced lung cancer [54]. GOT1 is mainly involved in amino acid metabolism and the urea and tricarboxylic acid cycle [55, 56]. Another study confirmed that, in OSCC, GOT1 is related to tumor invasion as well as shorter survival [57]. In addition, GOT1 is also closely related to esophageal squamous cell carcinoma [58] and can be used as a biomarker of prostate cancer [59]. However, the relationship among POLD2, HADHB, POLE3, ATIC, MGST1, GNPDA1 and OSCC remains unclear. POLD2 and POLE3 are necessary for DNA replication [60]. The expression of PLOD2 has been discovered to be increased in patients with bladder cancer and ovarian cancer, which predicts dismal patient survival [61, 62]. HADHB is related to fatty acid metabolism and has been reported to be elevated in renal clear cell carcinoma and colorectal cancer [63, 64]. The expression of ATIC is abnormally upregulated in hepatocellular carcinoma patients, which has a potential tumor-promoting effect [65]. Overexpression of MGST1 in tumor tissue can inhibit tumor cell apoptosis by inhibiting apoptosis-related signaling [66]. GNPDA1 is mainly related to glycolysis and amino acid metabolism; studies have shown that GNPDA1 is overexpressed in hepatocellular carcinoma [67]. It is speculated that the above genes could become prognostic markers of OSCC.
According to the results of these studies, we hypothesized that our identified gene signature composed of 11 MRGs could accurately predict OSCC prognosis. According to KM curve for training cohort, high-risk patients had markedly decreased OS compared with low-risk patients. The AUC of ROC curve for the risk score were 0.63, 0.70, 0.76 at 1, 3, and 5 years. A prediction models constructed by four autophagy-related genes perform well in predicting the overall survival of OSCC. The AUC values of these four genes was 0.65 [68]. Another study build a prognostic model for OSCC using 7 genes related to tumor mutational burden. From time-dependent ROC analysis, the AUC of 1, 3, and 5 years survivals of patients in the TCGA database were 0.67, 0.67 and 0.64 [69]. It indicated that this score predicted OSCC survival and that our as-constructed model performed well in predicting OSCC prognostic outcome. According to univariate as well as multivariate Cox regression on clinical features and risk scores of the patients, risk scores was the factor independently predicting prognosis. The 11 MRGs model performed well in prognosis prediction compared with additional clinicopathological characteristics alone. The above results were verified in the verification set, illustrating that our model has good universality and reliability. Finally, we established a nomogram by combining the risk model and the characteristics of clinical cases. The ROC curve and calibration curve both showed that our nomogram performed well in predicting OSCC prognosis. These findings suggest that combining risk score and additional clinical factors contributes to the accurate prediction of patient survival. As suggested by GSEA, the genes related to the regulation of metabolism were more highly associated with high-risk patients, which indicated the potent impact and regulation of MRGs on high-risk patients compared with low-risk counterparts. We also evaluated the relationship between immune infiltrating cells and risk scores, and the results suggest that the prognosis of the low-risk group seems to be better than that of the high-risk group. In addition, the results of single-cell sequencing and cell type identification show that the genes that make up OSCC prognosis model are limited to tumor cells rather than other cells. Collectively, our as-constructed 11 MRGs-based prognosis model may be adopted for predicting OSCC survival and providing more individualized treatment to OSCC patients.
MRGs can predict the prognosis of OSCC to some extent and are helpful for guiding patient treatment decision-making as well as follow-up visits. However, our study has some limitations. The data used in this study were acquired from public databases, so we cannot guarantee the integrity of patient data. In addition, due to the retrospective design of the study, there are potential biases related to imbalanced clinical characteristics. More experimental research and prospective works should be conducted for verifying MRGs’ prognosis prediction ability in OSCC.
Conclusion
In this study, we screened 11 MRGs related to the prognosis of OSCC through a variety of bioinformatics methods. Based on these 11 MRGs, a nomogram combining age, gender, staging and risk score was established. This nomogram can help clinicians more accurately determine the prognosis of OSCC and provide more individualized treatment to these patients.
Availability of data and materials
The datasets used and/or analysed during the current study are available from NCBI Gene Expression Omnibus (GEO: GSE41613) https://www.ncbi.nlm.nih.gov/geo/ and the cancer genome database (TCGA: OSCC) https://portal.gdc.cancer.gov/.
Abbreviations
- OC:
-
Oral cancer
- OSCC:
-
Oral squamous cell carcinoma
- AJCC:
-
American Joint Commission on Cancer
- MRGs:
-
Metabolism-related genes
- TCGA:
-
The Cancer Genome Atlas
- GEO:
-
Gene Expression Omnibus
- OS:
-
Overall survival
- HRs:
-
Hazard ratios
- CIs:
-
Intervals
- KM:
-
Kaplan–Meier
- ROC:
-
Receiver operating characteristic
- AUC:
-
Area under the curve
- GO:
-
Gene Ontology
- KEGG:
-
Kyoto Gene and Genomic Encyclopedia
References
Law ZJ, Khoo XH, Lim PT, Goh BH, Ming LC, Lee WL, Goh HP. Extracellular vesicle-mediated chemoresistance in oral squamous cell carcinoma. Front Mol Biosci. 2021;8:629888.
Bugshan A, Farooq I. Oral squamous cell carcinoma: metastasis, potentially associated malignant disorders, etiology and recent advancements in diagnosis. F1000Res. 2020;9:229.
Panarese I, Aquino G, Ronchi A, Longo F, Montella M, Cozzolino I, Roccuzzo G, Colella G, Caraglia M, Franco R. Oral and Oropharyngeal squamous cell carcinoma: prognostic and predictive parameters in the etiopathogenetic route. Expert Rev Anticancer Ther. 2019;19(2):105–19.
Leemans CR, Snijders P, Brakenhoff RH. The molecular landscape of head and neck cancer. Nat Rev Cancer. 2018;18(5):269–82.
Kondoh N, Mizuno-Kamiya M, Umemura N, Takayama E, Kawaki H, Mitsudo K, Muramatsu Y, Sumitomo S. Immunomodulatory aspects in the progression and treatment of oral malignancy. Jpn Dent Sci Rev. 2019;55(1):113–20.
Huang WC, Jang TH, Tung SL, Yen TC, Chan SH, Wang LH. A novel miR-365-3p/EHF/keratin 16 axis promotes oral squamous cell carcinoma metastasis, cancer stemness and drug resistance via enhancing beta5-integrin/c-met signaling pathway. J Exp Clin Cancer Res. 2019;38(1):89.
Sasahira T, Kirita T. Hallmarks of cancer-related newly prognostic factors of oral squamous cell carcinoma. INT J MOL SCI. 2018;19(8):7400.
Agarbati S, Mascitti M, Paolucci E, Togni L, Santarelli A, Rubini C, Fazioli F. Prognostic relevance of macrophage phenotypes in high-grade oral tongue squamous cell carcinomas. Appl Immunohistochem Mol Morphol. 2021;29(5):359–65.
Russo D, Mariani P, Caponio V, Lo RL, Fiorillo L, Zhurakivska K, Lo ML, Laino L, Troiano G. Development and validation of prognostic models for oral squamous cell carcinoma: a systematic review and appraisal of the literature. Cancers (Basel). 2021;13(22):996.
Moeckelmann N, Ebrahimi A, Tou YK, Gupta R, Low TH, Ashford B, Ch’Ng S, Palme CE, Clark JR. Prognostic implications of the 8th edition American Joint Committee on Cancer (AJCC) staging system in oral cavity squamous cell carcinoma. Oral Oncol. 2018;85:82–6.
Fu Z, Yu B, Liu M, Wu B, Hou Y, Wang H, Jiang Y, Zhu D. Construction of a prognostic signature in Ewing’s sarcoma: based on metabolism-related genes. Transl Oncol. 2021;14(12):101225.
Wang Z, Jiang Q, Dong C. Metabolic reprogramming in triple-negative breast cancer. Cancer Biol Med. 2020;17(1):44–59.
Ahmad F, Cherukuri MK, Choyke PL. Metabolic reprogramming in prostate cancer. Br J Cancer. 2021;125(9):1185–96.
Xu Z, Yuan KF. Lipid metabolic reprogramming and metabolic stress in liver cancer. Sichuan Da Xue Xue Bao Yi Xue Ban. 2021;52(4):561–5.
Chakraborty S, Balan M, Sabarwal A, Choueiri TK, Pal S. Metabolic reprogramming in renal cancer: events of a metabolic disease. Biochim Biophys Acta Rev Cancer. 2021;1876(1):188559.
Sun L, Suo C, Li ST, Zhang H, Gao P. Metabolic reprogramming for cancer cells and their microenvironment: beyond the Warburg Effect. Biochim Biophys Acta Rev Cancer. 2018;1870(1):51–66.
Luo T, Li Y, Nie R, Liang C, Liu Z, Xue Z, Chen G, Jiang K, Liu ZX, Lin H, et al. Development and validation of metabolism-related gene signature in prognostic prediction of gastric cancer. Comput Struct Biotechnol J. 2020;18:3217–29.
Kobayashi Y, Banno K, Kunitomi H, Takahashi T, Takeda T, Nakamura K, Tsuji K, Tominaga E, Aoki D. Warburg effect in Gynecologic cancers. J Obstet Gynaecol Res. 2019;45(3):542–8.
Furuta E, Okuda H, Kobayashi A, Watabe K. Metabolic genes in cancer: their roles in tumor progression and clinical implications. Biochim Biophys Acta. 2010;1805(2):141–52.
Jia D, Lu M, Jung KH, Park JH, Yu L, Onuchic JN, Kaipparettu BA, Levine H. Elucidating cancer metabolic plasticity by coupling gene regulation with metabolic pathways. Proc Natl Acad Sci U S A. 2019;116(9):3909–18.
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.
Long W, Gong X, Yang Y, Yang K. Downregulation of PER2 promotes tumor progression by enhancing glycolysis via the phosphatidylinositol 3-kinase/protein kinase b pathway in oral squamous cell carcinoma. J Oral Maxillofac Surg. 2020;78(10):1780–1.
Yang X, Wang Z, Li X, Liu B, Liu M, Liu L, Chen S, Ren M, Wang Y, Yu M, et al. SHMT2 desuccinylation by SIRT5 drives cancer cell proliferation. Cancer Res. 2018;78(2):372–86.
Chen G, Liu H, Zhang Y, Liang J, Zhu Y, Zhang M, Yu D, Wang C, Hou J. Silencing PFKP inhibits starvation-induced autophagy, glycolysis, and epithelial mesenchymal transition in oral squamous cell carcinoma. Exp Cell Res. 2018;370(1):46–57.
Adamski J, Suhre K. Metabolomics platforms for genome wide association studies–linking the genome to the metabolome. Curr Opin Biotechnol. 2013;24(1):39–47.
Huang J, Luo F, Shi M, Luo J, Ma C, Li S, Wei Y, Guo R, Li T. Construction and validation of a metabolic gene-associated prognostic model for cervical carcinoma and the role on tumor microenvironment and immunity. Aging (Albany NY). 2021;13(23):25072–88.
Wu Y, Wei X, Feng H, Hu B, Liu B, Luan Y, Ruan Y, Liu X, Liu Z, Wang S, et al. An eleven metabolic gene signature-based prognostic model for clear cell renal cell carcinoma. Aging (Albany NY). 2020;12(22):23165–86.
Cui L, Xue H, Wen Z, Lu Z, Liu Y, Zhang Y. Prognostic roles of metabolic reprogramming-associated genes in patients with hepatocellular carcinoma. Aging (Albany NY). 2020;12(21):22199–219.
Wu X, Yao Y, Li Z, Ge H, Wang D, Wang Y. Identification of a transcriptional prognostic signature from five metabolic pathways in oral squamous cell carcinoma. FRONT ONCOL. 2020;10:572919.
Verneau J, Sautes-Fridman C, Sun CM. Dendritic cells in the tumor microenvironment: prognostic and theranostic impact. Semin Immunol. 2020;48:101410.
Aponte-Lopez A, Munoz-Cruz S. Mast cells in the tumor microenvironment. Adv Exp Med Biol. 2020;1273:159–73.
Saikishore R, Velmurugan P, Ranjithkumar D, Latha R, Sathiamoorthi T, Arun A, Ravi AV, Sivakumar S. The circular RNA-miRNA axis: a special RNA signature regulatory transcriptome as a potential biomarker for OSCC. Mol Ther Nucleic Acids. 2020;22:352–61.
Yang Z, Liang X, Fu Y, Liu Y, Zheng L, Liu F, Li T, Yin X, Qiao X, Xu X. Identification of AUNIP as a candidate diagnostic and prognostic biomarker for oral squamous cell carcinoma. EBioMedicine. 2019;47:44–57.
Tian Z, Wang Z, Chen Y, Qu S, Liu C, Chen F, Ma L, Zhu J. Bioinformatics analysis of prognostic tumor microenvironment-related genes in the tumor microenvironment of hepatocellular carcinoma. Med Sci Monit. 2020;26:e922159.
Seborova K, Vaclavikova R, Soucek P, Elsnerova K, Bartakova A, Cernaj P, Bouda J, Rob L, Hruda M, Dvorak P. Association of ABC gene profiles with time to progression and resistance in ovarian cancer revealed by bioinformatics analyses. Cancer Med. 2019;8(2):606–16.
Hou C, Cai H, Zhu Y, Huang S, Song F, Hou J. Development and validation of autophagy-related gene signature and nomogram for predicting survival in oral squamous cell carcinoma. Front Oncol. 2020;10:558596.
Huang GZ, Lu ZY, Rao Y, Gao H, Lv XZ. Screening and identification of autophagy-related biomarkers for oral squamous cell carcinoma (OSCC) via integrated bioinformatics analysis. J Cell Mol Med. 2021;25(9):4444–54.
Li H, Zhang X, Yi C, He Y, Chen X, Zhao W, Yu D. Ferroptosis-related gene signature predicts the prognosis in Oral squamous cell carcinoma patients. BMC Cancer. 2021;21(1):835.
Zhang Z, Gao Z, Rajthala S, Sapkota D, Dongre H, Parajuli H, Suliman S, Das R, Li L, Bindoff LA, et al. Metabolic reprogramming of normal oral fibroblasts correlated with increased glycolytic metabolism of oral squamous cell carcinoma and precedes their activation into carcinoma associated fibroblasts. Cell Mol Life Sci. 2020;77(6):1115–33.
Ueda S, Goto M, Hashimoto K, Hasegawa S, Imazawa M, Takahashi M, Oh-Iwa I, Shimozato K, Nagao T, Nomoto S. Salivary CCL20 level as a biomarker for oral squamous cell carcinoma. Cancer Genom Proteom. 2021;18(2):103–12.
de Sa AM, de Sa RN, Bandeira CM, Chagas J, Pascoal M, Nepomuceno G, Da SMH, Alves M, Mendes MA, Dias M, et al. Identification of possible salivary metabolic biomarkers and altered metabolic pathways in South American patients diagnosed with oral squamous cell carcinoma. Metabolites. 2021;11(10):668.
Zeng Y, Zhang J, Xu M, Chen F, Zi R, Yue J, Zhang Y, Chen N, Chin YE. Roles of mitochondrial serine hydroxymethyltransferase 2 (SHMT2) in human carcinogenesis. J Cancer. 2021;12(19):5888–94.
Wu ZZ, Wang S, Yang QC, Wang XL, Yang LL, Liu B, Sun ZJ. Increased expression of SHMT2 is associated with poor prognosis and advanced pathological grade in oral squamous cell carcinoma. Front Oncol. 2020;10:588530.
Liao Y, Wang F, Zhang Y, Cai H, Song F, Hou J. Silencing SHMT2 inhibits the progression of tongue squamous cell carcinoma through cell cycle regulation. Cancer Cell Int. 2021;21(1):220.
Zhang P, Yang Q. Overexpression of SHMT2 predicts a poor prognosis and promotes tumor cell growth in bladder cancer. Front Genet. 2021;12:682856.
Jin M, Lee WK, You MH, Jang A, Cheng SY, Kim WG, Jeon MJ, Lee YM. SHMT2 expression as a diagnostic and prognostic marker for thyroid cancer. Endocr Connect. 2021;10(6):630–6.
de Campos RP, Schultz IC, de Andrade MP, Davies S, Gasparin MS, Bertoni A, Buffon A, Wink MR. Cervical cancer stem-like cells: systematic review and identification of reference genes for gene expression. Cell Biol Int. 2018;42(2):139–52.
Wu T, Jiao Z, Li Y, Su X, Yao F, Peng J, Chen W, Yang A. HPRT1 promotes chemoresistance in oral squamous cell carcinoma via activating MMP1/PI3K/Akt signaling pathway. Cancers (Basel). 2022;14(4):668.
Li S, Mai Z, Gu W, Ogbuehi AC, Acharya A, Pelekos G, Ning W, Liu X, Deng Y, Li H, et al. Molecular subtypes of oral squamous cell carcinoma based on immunosuppression genes using a deep learning approach. Front Cell Dev Biol. 2021;9:687245.
Zhang M, Zeng X, Yang Q, Xu J, Liu Z, Zhou Y, Cao Y, Zhang X, An X, Xu Y, et al. Ablation of myeloid ADK (adenosine kinase) epigenetically suppresses atherosclerosis in ApoE(−/−) (Apolipoprotein E Deficient) mice. Arterioscler Thromb Vasc Biol. 2018;38(12):2780–92.
Shamloo B, Kumar N, Owen RH, Reemmer J, Ost J, Perkins RS, Shen HY. Dysregulation of adenosine kinase isoforms in breast cancer. Oncotarget. 2019;10(68):7238–50.
Rai B, Kaur J, Jacobs R, Anand SC. Adenosine deaminase in saliva as a diagnostic marker of squamous cell carcinoma of tongue. Clin Oral Investig. 2011;15(3):347–9.
Kelgandre DC, Pathak J, Patel S, Ingale P, Swain N. Adenosine deaminase-a novel diagnostic and prognostic biomarker for oral squamous cell carcinoma. Asian Pac J Cancer Prev. 2016;17(4):1865–8.
Zanini D, Manfredi LH, Pelinson LP, Pimentel VC, Cardoso AM, Carmo AGV, Santos C, Gutierres JM, Morsch VM, Leal D, et al. ADA activity is decreased in lymphocytes from patients with advanced stage of lung cancer. Med Oncol. 2019;36(9):78.
Kremer DM, Nelson BS, Lin L, Yarosz EL, Halbrook CJ, Kerk SA, Sajjakulnukit P, Myers A, Thurston G, Hou SW, et al. GOT1 inhibition promotes pancreatic cancer cell death by ferroptosis. Nat Commun. 2021;12(1):4860.
Wang Q, Zhang Q, Luan S, Yang K, Zheng M, Li K, Chen L, Li H. Adapalene inhibits ovarian cancer ES-2 cells growth by targeting glutamic-oxaloacetic transaminase 1. Bioorg Chem. 2019;93:103315.
Busso-Lopes AF, Carnielli CM, Winck FV, Patroni F, Oliveira AK, Granato DCECR, Domingues RR, Pauletti BA, Riano-Pachon DM, et al. A reductionist approach using primary and metastatic cell-derived extracellular vesicles reveals hub proteins associated with oral cancer prognosis. Mol Cell Proteom. 2021;20:100118.
Zhou S, Guo Z, Lv X, Zhang X. CircGOT1 promotes cell proliferation, mobility, and glycolysis-mediated cisplatin resistance via inhibiting its host gene GOT1 in esophageal squamous cell cancer. Cell Cycle. 2022;21(3):247–60.
Lee B, Mahmud I, Marchica J, Derezinski P, Qi F, Wang F, Joshi P, Valerio F, Rivera I, Patel V, et al. Integrated RNA and metabolite profiling of urine liquid biopsies for prostate cancer biomarker discovery. Sci Rep. 2020;10(1):3716.
Bellelli R, Belan O, Pye VE, Clement C, Maslen SL, Skehel JM, Cherepanov P, Almouzni G, Boulton SJ. POLE3-POLE4 Is a Histone H3–H4 chaperone that maintains chromatin integrity during DNA replication. Mol Cell. 2018;72(1):112–26.
Givechian KB, Garner C, Garban H, Rabizadeh S, Soon-Shiong P. CAD/POLD2 gene expression is associated with poor overall survival and chemoresistance in bladder urothelial carcinoma. Oncotarget. 2018;9(51):29743–52.
Zhu Y, Lu H, Zhang D, Li M, Sun X, Wan L, Yu D, Tian Y, Jin H, Lin A, et al. Integrated analyses of multi-omics reveal global patterns of methylation and hydroxymethylation and screen the tumor suppressive roles of HADHB in colorectal cancer. Clin Epigenetics. 2018;10:30.
Zhao Z, Liu Y, Liu Q, Wu F, Liu X, Qu H, Yuan Y, Ge J, Xu Y, Wang H. The mRNA expression signature and prognostic analysis of multiple fatty acid metabolic enzymes in clear cell renal cell carcinoma. J Cancer. 2019;10(26):6599–607.
Li M, Jin C, Xu M, Zhou L, Li D, Yin Y. Bifunctional enzyme ATIC promotes propagation of hepatocellular carcinoma by regulating AMPK-mTOR-S6 K1 signaling. Cell Commun Signal. 2017;15(1):52.
Kuang F, Liu J, Xie Y, Tang D, Kang R. MGST1 is a redox-sensitive repressor of ferroptosis in pancreatic cancer cells. Cell Chem Biol. 2021;28(6):765–75.
Torti SV, Torti FM. Iron: the cancer connection. Mol Aspects Med. 2020;75:100860.
Li D, Cheng X, Zheng W, Chen J. Glucosamine-6-phosphate isomerase 1 promotes tumor progression and indicates poor prognosis in hepatocellular carcinoma. Cancer Manag Res. 2020;12:4923–35.
Zhu L, Yan D, Chen Y, Chen S, Chen N, Han J. The identification of autophagy-related genes in the prognosis of oral squamous cell carcinoma. Oral Dis. 2020;26(8):1659–67.
Wu F, Du Y, Hou X, Cheng W. A prognostic model for oral squamous cell carcinoma using 7 genes related to tumor mutational burden. BMC Oral Health. 2022;22(1):152.
Acknowledgements
We thank Dr.Jianming Zeng (University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes. We thank Xiaoqi Wu (Genergy Biotechnology Shanghai Co., Ltd) for his help in bioinformatics analysis.
Funding
This research was supported by the Science and Technology Foundation of Shandong Province (ZR2021MH033), Traditional Chinese Medicine Science and Technology Project of Shandong Province (2021M024) and China Postdoctoral Science Foundation (No. 2018M632679).
Author information
Authors and Affiliations
Contributions
LL and ZC conceived and designed the study. JZ, CM, HQ, XH and JL collected and analyzed the data, and wrote the paper. ZW, CZ and XL downloaded and organized data, and drew the figures. All authors read and approved the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
All 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.
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.
About this article
Cite this article
Zhang, J., Ma, C., Qin, H. et al. Construction and validation of a metabolic-related genes prognostic model for oral squamous cell carcinoma based on bioinformatics. BMC Med Genomics 15, 269 (2022). https://doi.org/10.1186/s12920-022-01417-3
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12920-022-01417-3