A Novel Autophagy-Related Prognostic Risk Model and a Nomogram for Survival Prediction of Oral Cancer Patients
Background. This study is aimed at constructing a risk signature to predict survival outcomes of ORCA patients. Methods. We identified differentially expressed autophagy-related genes (DEARGs) based on the RNA sequencing data in the TCGA database; then, four independent survival-related ARGs were identified to construct an autophagy-associated signature for survival prediction of ORCA patients. The validity and robustness of the prognostic model were validated by clinicopathological data and survival data. Subsequently, four independent prognostic DEARGs that composed the model were evaluated individually. Results. The expressions of 232 autophagy-related genes (ARGs) in 127 ORCA and 13 control tissues were compared, and 36 DEARGs were filtered out. We performed functional enrichment analysis and constructed protein–protein interaction network for 36 DEARGs. Univariate and multivariate Cox regression analyses were adopted for searching prognostic ARGs, and an autophagy-associated signature for ORCA patients was constructed. Eventually, 4 desirable independent survival-related ARGs (WDR45, MAPK9, VEGFA, and ATIC) were confirmed and comprised the prognostic model. We made use of multiple ways to verify the accuracy of the novel autophagy-related signature for survival evaluation, such as receiver-operator characteristic curve, Kaplan–Meier plotter, and clinicopathological correlational analyses. Four independent prognostic DEARGs that formed the model were also associated with the prognosis of ORCA patients. Conclusions. The autophagy-related risk model can evaluate OS for ORCA patients independently since it is accurate and stable. Four prognostic ARGs that composed the model can be studied deeply for target treatment.
Death rates are increasing for oropharyngeal cancer overall by 0.5% per year from 2009 to 2018. There is an estimated number of 54,010 new cases of oropharyngeal cancer (accounting for about 2.8% of new cancers), and 10,850 will die of it (accounting for about 1.8% of all deaths from cancer) in 2021 in the United States . Among all oropharyngeal cancer patients, approximately half of those occur specifically in the oral cavity and were called oral cancer (ORCA) . All over the world, 377,713 people were diagnosed with oral cancer in 2020, and the cancer caused 177,757 deaths; it means there was one person who died due to ORCA every 3 minutes . The treatment and prognosis of ORCA patients relied on the staging system and conventional prognostic factors used in clinical practice, for example, the tumor-node-metastasis (TNM) staging system was the most important and well-known prognostic factor for ORCA patients , but the staging system was not accurate enough and prognostic factors were inadequate and nonspecific, leading to discouraging overall prognosis. Recent statistics show that the overall 5-year survival rate for ORCA patients was 52.0% in China . Therefore, searching molecular prognostic markers which could be incorporated into the prognostic system is significant; they can benefit predicting clinical outcomes and targeted therapy. This study attempts to refine the prognostic signature of ORCA and used it in the clinical setting.
Autophagy is a process that can circulation and reuse cellular components; it is nonspecific. Autophagy-related genes (ARGs) are molecules that participate in autophagy. In recent studies, many researchers deemed that autophagy is an important process in ORCA [6–8]. Most of the molecular mechanisms which are involved in autophagy regulation deeply participate in signaling pathways that assumed significant roles in ORCA control [9, 10]. For example, autophagy is a significant AKT/mTOR pathway target, and the antioncogene that negatively regulates mTOR (such as AMPK and PTEN) will activate autophagy. In addition, many autophagy-inducing proteins are either oncoproteins or tumor suppressor proteins . Some studies have proven that the autophagy inhibitor can enhance efficacy against aggressive ORCA [12, 13]. However, how to predict the prognosis of ORCA patients by several ARGs is still unclear. Therefore, our study utilized several screened ARGs to establish the prognostic risk model of ORCA. In this study, the relationship between prognostic ARGs and clinicopathological features in 127 ORCA patients were evaluated, and an autophagy-related signature was constructed to predict survival for ORCA patients. The risk score model had been verified from several aspects, and we proved its accuracy; we also verified the prognosis value of each ARG that comprised the autophagy-related risk model and proved them as prognostic biomarkers in ORCA; our report may shed light on evaluating the prognosis and targeting treatment of ORCA.
2.1. Data Acquisition
In total, 232 genes presently known as autophagy-related genes (ARGs) were downloaded from HADb (Human Autophagy Database, http://autophagy.lu/). RNA-seq data and clinical features of 127 ORCA patients and 13 nontumor samples were obtained from the TCGA website (The Cancer Genome Atlas database, https://tcga-data.nci.nih.gov/tcga/).
2.2. Identification of DEARGs and Functional Annotation
The “limma” package in R was applied to compare differential expressions between ORCA and nontumor counterparts; the criteria is false discovery rate and . Then, we performed enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for DEARGs to find the major biological attributes. The R packages of ClusterProfiler and org.Hs.eg.db were used to perform GO and KEGG function enrichment analyses. The visualization of results was implemented with the R “ggplot2” package. STRING database (https://string-db.org/) was used to construct the protein–protein interaction (PPI) network of DEARGs. Cytoscape software visualized the interaction of the PPI network.
2.3. Construction and Verification of ARG-Based Prognostic Risk Model
Univariate Cox proportional risk regression analysis was performed to screen prognostic ARGs. All statistically significant prognostic ARGs were selected as candidates for multivariate Cox regression analysis and used to construct the prognostic risk model. The multivariate analysis provided the regression coefficient of prognostic ARGs in the risk score model.
The median risk score was determined as the threshold to divide the ORCA patients into the high- or low-risk group. Four risk genes comprising the risk model were used to construct the nomogram. We constructed a nomogram based on 4 independent prognosis-related ARGs to predict survival rates of patients at 1, 2, and 3 years via the rms package in R 4.0.1.
The predictive accuracy of the risk models was evaluated by the receiver-operator characteristic (ROC) analysis; besides, we compared survival probability of high-risk and low-risk groups using the “ggplot2” package in R. We also assessed the association between the prognostic risk model and classical clinical parameters (such as tumor grade, pathological stage, and T/N classification) by Cox proportional hazard regression analysis.
2.4. Validation for Prognostic Value of 4 Risk ARGs which Formed the Model
The relationships between expressions of 4 independent prognostic-related ARGs (WDR45, MAPK9, VEGFA, and ATIC) and OS of ORCA patients were evaluated individually by Kaplan–Meier survival curves. We also observed the mRNA expression pattern of the 4 genes in ORCA and normal tongue tissues via the Oncomine database (https://www.oncomine.org/).
To explore the roles of 4 risk ARGs in ORCA, a single-gene gene set enrichment analysis (GSEA) was performed. GSEA was conducted using a molecular signature database’s (MSigDB’s) c2.cp.kegg.v7.0.symbols.gmt gene sets in GSEA 3.0 software to identify gene sets that were significantly correlated with the expression of risk genes. We exhibited the top 10 biological processes most associated with the expression of the 4 risk genes.
3.1. Identification of DEARGs
The expression profile of 232 ARGs in 127 ORCA tissue samples and 13 normal mouth samples were compared and analyzed. Finally, there were 36 differentially expressed ARGs including 29 upregulated ARGs and 7 downregulated ARGs with a threshold of . The volcano plot, heat map, and box plot visualized the expression pattern of all DEARGs between ORCA and normal tissues (Figures 1(a)–1(c)). We listed all DEARGs in Table 1, including the log2FoldChange and statistical significance.
3.2. Functional Annotation for All DEARGs
In Figure2(a), all DEARGs are linked to form a protein-protein interaction network. The DEARGs are enriched in several biological processes (BP) such as autophagy and protein localization. In the molecular function (MF) term of GO analysis, some protein binding-related functions including “protease binding,” “integrin binding,” and “ubiquitin protein ligase binding” were significant roles these DEARGs played. For cellular components (CC), “focal adhesion,” “vacuolar membrane,” and “cell substrate” were significantly enriched (Figure2(b)). KEGG analysis exhibited that the DEARGs were significantly associated with apoptosis, platinum drug resistance, and hepatitis C (Figure2(c)).
3.3. Establishment of Prognosis Prediction Model with ARGs
We performed the univariate Cox proportional hazard analysis on all ARGs; a total of 13 ARGs were significantly associated with the prognosis of ORCA patients () (Figure 3(a)). Then, these prognostic ARGs entered the multivariate Cox regression analysis, and finally, 4 of them (WDR45, MAPK9, VEGFA, and ATIC) were screened as independent prognostic genes to construct the prognosis prediction model (Figure 3(b)). For each ORCA patient, the . We divided the 126 ORCA cases into high- or low-risk groups according to the median values of the risk score. In Figure 3(c), the above 4 independent prognostic ARGs were incorporated into a nomogram model for predicting the individualized probability of survival times in clinical practice in ORCA patients. The score of every ARG can be found through the point scale located at the top of the nomogram. Then, the points of each ARG were summed, thereby estimating survival probability at 1, 2, and 3 years.
3.4. Testing the Risk Model
Patients with ORCA were stratified into high- and low-risk groups based on the median risk score, and Kaplan-Meier plots indicated that patients in the low-risk group had a significantly better prognosis () (Figure 4(a)). The results of ROC curves revealed that this prognostic risk model had good predictive performance in predicting the overall survival of ORCA patients with areas under the ROC curve (AUCs) of 0.742 (Figure 4(b)). Figures 4(c)–4(e) indicate the risk distribution of patients, survival time of patients in high- and low-risk groups, and expression profile heat map of the 4 ARGs that formed the risk model.
We grouped the patients by age, sex, histological grade, clinical stage, and T/N classification to explore the association between the risk model and clinicopathological features. In Figure 5, we can find a trend of risk scores which were higher in patients classified with high grade, advanced pathologic stage, and terminal T/N stage (Figure 5(a)). The expressions of MAPK9, VEGFA, and ATIC all had a trend of being upregulated in patients with high grade, advanced pathologic stage, and terminal T/N stage; in addition, the expressions of WDR45 had an opposite trend (Figures 5(b)–5(e)). The expression regularity of 4 independent prognostic ARGs that comprised the risk model in patients with different clinicopathological features was in accordance with their expression trend in different risk groups. The expressions of MAPK9, VEGFA, and ATIC were higher in the high-risk group, and the expressions of WDR45 was higher in the low-risk group (Figure 4(e)). It proved that the risk prognostic model was better than any other clinicopathological features for predicting the prognosis for ORCA patients.
The results in Table 2 suggested that our prognostic risk model could independently predict the prognosis of ORCA. Although T classification was also significantly associated with survival, ROC analysis demonstrated that T classification () was not as reliable as the risk score ().
3.5. Validation of 4 Risk Genes that Comprised the Risk Model
WDR45, MAPK9, VEGFA, and ATIC were 4 prognostic-related ARGs which formed the prognostic risk model. As shown in Figure 6(a), upregulation of WDR45 was strongly correlated with longer patient survival, and MAPK9, VEGFA, and ATIC overexpression reduced OS in patients with ORCA. We also verified the expression of 4 risk genes in the Oncomine database (Figure 6(b)) and found that the expression trend of 4 prognostic ARGs is in accordance with their expression patterns in ORCA, and there are no tumor tissues in previous results in this paper in the TCGA database as shown in Figure 4(e) and Figures 5(b)–5(e). Single-gene GSEA of the 4 risk genes revealed their potential function in ORCA (Figure 7).
Autophagy is a conserved and dynamic process whose function maintains cellular homeostasis . Multiple studies had confirmed that autophagy played a significant role in many cancers [15–17]. Important experimental evidences had proven autophagy’s potential as a therapeutic target for ORCA [18–20]. Hence, we constructed an autophagy-related prognostic risk model in this paper for prognosis prediction for ORCA patients.
The autophagy-related prognostic model is formed by 4 ARGs, including WDR45, MAPK9, VEGFA, and ATIC. According to our risk model, the expressions of MAPK9, VEGFA, and ATIC in the high-risk group were higher than those in the low risk-group; it is worth noting that their expressions in ORCA tissues were higher than those in the control, and their expressions all had a trend of upregulation in patients with advanced pathologic stage and high grade; in addition, high expressions of the 3 genes were significantly related to worse OS. What this means is that MAPK9, VEGFA, and ATIC were indeed high-risk factors in ORCA. The expression law of WDR45 was opposite to MAPK9, VEGFA, and ATIC and proved it was a low-risk factor in ORCA. The risk score model is significantly associated with clinicopathological indicators of ORCA patients. Our paper proved that the risk model comprised of WDR45, MAPK9, VEGFA, and ATIC for evaluating the prognosis of ORCA patients is clinically practicable.
MAPK9 is also called JNK2; recent studies have highlighted the oncogenic potential of MAPK9 in several human cancer cells, such as lung  and glioblastoma . The role of MAPK9 in oral cancer is still controversial . VEGFA encoded vascular endothelial growth factor A; many researchers reported that it was upregulated in oral squamous cell carcinomas [24, 25]. ATIC was considered an effective target for chemoradiosensitization . Although there are some researches about independent prognostic ARGs in recent years, they are not incorporated into a prognostic evaluation system in the right proportions and ways, so exploring an effective prognostic risk model and novel targets for ORCA therapy is necessary.
Functional enrichment analysis showed that 36 DEARGs were mainly involved in apoptosis and inflammation-associated pathways. Dwivedi et al. reported that to regulate apoptosis may be the best way against ORCA . A lot of reports deemed that oral inflammation promotes ORCA [28, 29]. Hence, clinical doctors maybe should pay more attention to the effect of inflammation on ORCA patients.
In this paper, we demonstrated that the novel prognostic risk model had excellent potential as a prognostic predictor which is even better than traditional clinicopathological indicators in ORCA. It can be incorporated into the clinical evaluation indicators to better predict clinical outcomes. Individual assessment of 4 risk genes that comprised the risk model in ORCA further proved that WDR45, MAPK9, VEGFA, and ATIC all play important parts in ORCA.
A novel autophagy-related model based on the expression levels of 4 risk genes was explored for survival prediction for ORCA patients. The 4 risk genes can benefit the underlying molecular mechanisms of ORCA and be utilized as potential therapeutic targets.
|ARGs:||Autophagy- related genes|
|AUC:||Area under the curve|
|DEARGs:||Differentially expressed autophagy- related genes|
|DEGs:||Differentially expressed genes|
|FDR:||False discovery rate|
|GSEA:||Gene set enrichment analysis|
|KEGG:||Kyoto Encyclopedia of Genes and Genomes|
|KM plotter:||Kaplan-Meier plotter|
|Log2FC:||Log2FoldChange/logarithm of fold change|
|TCGA:||The Cancer Genome Atlas.|
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.
Conflicts of Interest
The authors have no conflicts of interest to declare.
HF designed this study and drafted this paper. HF and XC conducted experiments and did statistical analysis. HF searched the references and participated in making the figure. HF revised the manuscript and supervised all the work.
The study was supported by the Nosocomial Scientific Research Fund Projects from the International Peace Maternity and Child Health Hospital of Shanghai Jiao Tong University School of Medicine (No. GFY5801), Clinical Research Special Projects from the Shanghai Municipal Health Commission (No. 20204Y0230), and Shanghai Sailing Program from the Shanghai Science and Technology Committee (No. 19YF1452200). We also appreciate the support of the Youth Science and Technology Innovation Studio, Shanghai Jiao Tong University School of Medicine.
Y. Wang, C. Wang, H. Tang et al., “Decrease of autophagy activity promotes malignant progression of tongue squamous cell carcinoma,” Journal of Oral Pathology & Medicine: official publication of the International Association of Oral Pathologists and the American Academy of Oral Pathology., vol. 42, no. 7, pp. 557–564, 2013.View at: Publisher Site | Google Scholar
S. Patil, R. S. Rao, and A. T. Raj, “Dual role of autophagy in oral cancer,” Journal of International Oral Health: JIOH, vol. 7, no. 6, pp. i–ii, 2015.View at: Google Scholar
H. Fei, S. Chen, and C. Xu, “Construction autophagy-related prognostic risk signature to facilitate survival prediction, individual treatment and biomarker excavation of epithelial ovarian cancer patients,” Journal of Ovarian Research, vol. 14, no. 1, p. 41, 2021.View at: Publisher Site | Google Scholar
C. B. Sales, M. E. Buim, R. O. de Souza et al., “Elevated VEGFA mRNA levels in oral squamous cell carcinomas and tumor margins: a preliminary study,” Journal of Oral Pathology & Medicine: official publication of the International Association of Oral Pathologists and the American Academy of Oral Pathology, vol. 45, no. 7, pp. 481–485, 2016.View at: Publisher Site | Google Scholar
T. Naruse, G. Kawasaki, S. Yanamoto, A. Mizuno, and M. Umeda, “Immunohistochemical study of VEGF expression in oral squamous cell carcinomas: correlation with the mTOR-HIF-1alpha pathway,” Anticancer Research, vol. 31, no. 12, pp. 4429–4437, 2011.View at: Google Scholar