An autophagy-related model of 4 key genes for predicting prognosis of patients with laryngeal cancer.

Autophagy, a major cause of cancer-related death, is correlated with the pathogenesis of various diseases including cancers. Our study aimed to develop an autophagy-related model for predicting prognosis of patients with laryngeal cancer.We analyzed the correlation between expression profiles of autophagy-related genes (ARGs) and clinical outcomes in 111 laryngeal cancer patients from The Cancer Genome Atlas (TCGA). Afterward, gene functional enrichment analyses of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed to find the major biological attributes. Univariate Cox regression analyses and multivariate Cox regression analyses were performed to screen ARGs whose expression profiles were significantly associated with laryngeal cancer patients overall survival (OS). Furthermore, to provide the doctors and patients with a quantitative method to perform an individualized survival prediction, we constructed a prognostic nomogram.Thirty eight differentially expressed ARGs were screened out in laryngeal cancer patients through the TCGA database. Related functional enrichments may act as tumor-suppressive roles in the tumorigenesis of laryngeal cancer. Subsequently, 4 key prognostic ARGs (IKBKB, ST13, TSC2, and MAP2K7) were identified from all ARGs by the Cox regression model, which significantly correlated with OS in laryngeal cancer. Furthermore, the risk score was constructed, which significantly divided laryngeal cancer patients into high- and low-risk groups. Integrated with clinical characteristics, gender, N and the risk score are very likely associated with patients OS. A prognostic nomogram of ARGs was constructed using the Cox regression model.Our study could provide a valuable prognostic model for predicting the prognosis of laryngeal cancer patients and a new understanding of autophagy in laryngeal cancer.

Introduction

Autophagy is a major catabolic system, of which cells can recycle organelles and proteins by degradation in the lysosomes.[ Autophagy selectively targets intracellular microbes, dysfunctional organelles, and pathogenic proteins, and deficiencies in these processes, which may lead to diseases, including cancers, neurodegenerative diseases, and inflammation.[ However, the mechanism of autophagy suppressing tumorigenesis is still unclear and inconclusive.[ Because of the complex function of autophagy in cancers, further researches may be focused on the relation of autophagy and tumors, especially biological processes.[ And then applying this knowledge to a well-designed therapeutic strategy could be a new method of cancer therapy.[ Even though plenty of researchers have made a lot of effort, the question of whether autophagy is a friend or a foe for cancers cannot be drawn a reliable conclusion for now.[

Laryngeal cancer, one of the most common head and neck cancers worldwide, has a poor prognosis and high mortality.[ The numbers of laryngeal cancer new cases and related deaths in the United States were estimated to be 12,410 and 3760 in 2019.[ Recently, many studies reported that autophagy inhibition could be a promising therapeutic target for laryngeal cancer, which provided a new route for the clinical management of laryngeal cancer.[ Therefore, exploring the related molecular biomarkers of autophagy would have an attractive value in estimating and treating laryngeal cancer, which may be an important therapeutic trend of laryngeal cancer.

In our study, we analyzed the correlation between expression profiles of autophagy-related genes (ARGs) and clinical outcomes in 111 laryngeal cancer patients from The Cancer Genome Atlas (TCGA). Afterward, gene functional enrichment analyses of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed to find the major biological attributes. Univariate Cox regression analyses and multivariate Cox regression analyses were performed to screen the ARGs whose expression profiles were significantly associated with laryngeal cancer patients overall survival (OS). To facilitate the utility of the selected prognostic ARGs in routine clinical practice, we integrated risk score with clinical factors to improve the prognostic efficiency of laryngeal cancer patients. By using the multivariate Cox regression analysis, the ARGs expressions based on age, gender, grade, stage, T, M, and N were analyzed to determine whether the ARG could be an independent predictor for the patients of laryngeal cancer or not. Furthermore, to provide the doctors and patients with a quantitative method to perform an individualized survival prediction, we constructed a prognostic nomogram. These findings could also provide an effective multi-dimensional biomarker strategy, which would be effective in monitoring autophagy and predicting the prognosis in laryngeal cancer patients.

Methods

Data acquisition

The Human Autophagy Database (HADb, http://www.autophagy.lu/index.html) is the first autophagy-dedicated database, which provides a complete and up-to-date list of human genes involved in autophagy.[ The list of ARGs was obtained from HADb. ARGs RNA-sequencing (RNA-seq) data and the clinical data of laryngeal cancer were downloaded and extracted from the TCGA database (https://portal.gdc.cancer.gov).[ The dataset contained 111 patients with solid tumor and 12 patients with non-tumor tissues. Because our study is a bioinformatics study, there is no need of ethics committee or institutional review board to approve the study.

Differentially expressed ARGs analysis

The “limma” package was applied in R statistical software to estimate differentially expressed ARGs between laryngeal cancer and non-tumor samples. A volcano plot and heat map were drawn using the “ggplot2” package and the “pheatmap” package of R software. Genes that showed at least 1-fold changes and an adjusted P value less than .05 were considered as the significant differentially expressed ARGs. The visualization of gene expression of the significant differentially expressed ARGs was drawn the “ggpubr” package of R software. Afterward, gene functional enrichment analyses of GO and KEGG were performed to find the major biological attributes.[ The packages of “clusterProfiler”, “org.Hs.eg.db”, “enrichplot”, “ggplot2”, and“GOplot”, were applied in R statistical software to conduct enriched GO and KEGG enrichment analyses and the visualization of enrichment terms.

Construction of risk score formula and ARGs-based prognostic model

Univariate Cox regression analyses were performed to screen ARGs whose expression profiles were significantly associated with laryngeal cancer patients OS. If the number of differentially expressed ARGs was low, we would use all the data of ARGs related to laryngeal cancer. Then, multivariate Cox regression analyses were conducted to remove the genes which might be considered not independent indicators. Subsequently, several prognostic ARGs were selected and the related risk score was developed.[ The formula of risk score is based on a linear combination of the relative expression level of genes multiplied regression coefficients, which represented the relative weight of genes. The risk score formula was built as following: Risk score = sum of coefficients × expression level of ARGs.[ By the risk cutoff value (the median risk score value), patients of laryngeal cancer were separated into high-risk and low-risk groups.[ The survival curve was plotted by Kaplan–Meier (K–M) method and the difference in the survival rates between high-risk and low-risk groups was assessed using the log-rank test.

By using the multivariate Cox regression analysis, ARGs expressions based on age, gender, grade, stage, T, M, and N were analyzed to determine whether the ARG could be an independent predictor of the patients with laryngeal cancer or not.

Construction of the nomogram and evaluation of ARGs-based prognostic model

To provide doctors and patients with a quantitative method to perform an individualized survival prediction, we constructed a prognostic nomogram of the ARGs using the Cox regression model.

To further assess the predictive performance of the ARGs-based prognostic model, we calculated the area under the curve (AUC). Furthermore, the predictive accuracy of the ARGs-based prognostic model was compared with other risk factors using receiver operating characteristic (ROC) analysis. Moreover, calibration curves were used to assess the agreements between model predicted outcomes and actual outcomes. Besides, the visualization of high-risk and low-risk groups was drawn using the “ggplot2” package of R software.

The validation of the autophagy-related model based on TCGA database

The autophagy-related model was validated using the TGGA database. Clinical data on survival and outcome were also downloaded from the TCGA database. The survival curves were plotted by K–M method.

Statistical analysis

All statistical analyses were performed using SPSS 24.0 (Chicago, IL, USA) and R version 3.5.3 (https://www.r-project.org/). R software was performed to draw plots. Univariate Cox regression analyses were calculated to evaluate the association between expression profiles and OS, and hazard ratios (HR) were calculated using the Cox regression model. The Multivariate Cox proportional hazards regression model was calculated to construct the ARGs-based prognostic model. ROC curve and the area under the ROC curve for each dataset were performed by the “survivalROC” package in R software. AUC > 0.7 was considered as a good performance. All tests were two-tailed and all statistical significance was defined as P < .05.

Results

Differential expression screening and identification of ARGs

Expression values of 222 ARGs were extracted. Based on the criteria of an FDR <0.05 and |log2(Fold Change)|>1, we finally obtained 29 up-regulated and 9 down-regulated ARGs (Fig. 1A, 1B). Moreover, visual scatter plots were exhibited the expression pattern of the 38 differentially expressed ARGs between laryngeal cancer and non-tumor tissues (Fig. 1C). Scatter plots showed expression patterns of ARGs.

Figure 1

Differentially expressed autophagy-related genes (ARGs) between laryngeal cancer and non-tumor tissues. (A) The volcano plot for the 222 ARGs from The Cancer Genome Atlas (TCGA) database. Red: upregulation; blue:downregulation. (B) Hierarchical clustering of differentially expressed ARGs expression levels. (C) The expression patterns of 38 ARGs.

Functional annotation of differentially expressed ARGs

The functional enrichment analysis of 38 differentially expressed ARGs was performed. GO terms function and KEGG pathway enrichment of these differentially expressed ARGs were summarized in Table 1. Using the packages of “clusterProfiler” and “org.Hs.eg.db” in R software, we proved that the top enriched GO terms for biological processes were: neuron apoptotic process, neuron death, and autophagy; for cellular components were: integrin complex, protein complex involved in cell adhesion, and lysosomal lumen; and for molecular function: fibronectin binding, protein tyrosine kinase activity, and receptor-ligand activity. The overview visualization of the GO analysis results is displayed in Figure 2A and B. Moreover, KEGG pathway enrichment analysis for the differentially expressed ARGs showed that they are notably associated with EGFR tyrosine kinase inhibitor resistance, apoptosis, ErbB signaling pathway, and so on. As shown in Figure 2C, the high Z-score of enriched pathways indicated that related pathways were more likely to be increased. The heatmap and circle plot of the relationship between ARGs and pathways were visually displayed (Fig. 2D and E).

Table 1

GO and KEGG analysis of differentially expressed autophagy-related genes.

Figure 2

Functional enrichment analyses of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). (A) The bubble plot of enriched GO terms. (B) The circos plot of enriched GO terms. (C) The circular scatter plot of enriched KEGG terms. (D) The heatmap of the relationship between ARGs and KEGG pathways. (E) The circos plot of TOP 10 enriched GO terms.

Identification of prognostic ARGs

The relationships between the expression profiles of ARGs and OS were assessed based on the data obtained from TCGA. Using univariate Cox regression analyses, 14 prognosis-related ARGs were selected with the criteria of a P < .05 (Table 2, Fig. 3A). To improve the robustness, the further multivariate Cox regression model was conducted (Table 2, Fig. 3B). Therefore, 4 genes including IKBKB, ST13, TSC2, and MAP2K7 were identified (Table 2). The results from K-M analysis indicated that the overexpression of ST13 was strongly correlated with the worse OS of laryngeal cancer patients (HR = 3.685, 95% CI = 1.571–8.645, P = .003; Fig. 3B). On the contrary, Figure 3B showed that up-regulated IKBKB, up-regulated TSC2, and up-regulated MAP2K7 indicated laryngeal cancer patients may likely have a longer survival time.

Table 2

Expression and Cox regression analysis data of the prognosis-related ARGs in laryngeal cancer by TCGA.

Figure 3

Forest plots of autophagy-related genes (ARGs) and autophagy-related risk score of laryngeal cancer patients. (A) Univariate analyses of ARGs. (B) Multivariate analyses of 4 key ARGs. (C) Kaplan–Meier plot represents that patients in the high-risk group had a significantly worse overall survival (OS) than those in the low-risk group. (D) The risk score distribution of patients in the TCGA dataset. (E) The number of patients in different risk groups. (F) The survival time of patients in the TCGA dataset. (G) The heatmap of the 4 key ARGs expression profiles in the TCGA dataset.

Development of risk score formula and ARGs-based prognostic model

To facilitate the utility of the selected prognostic ARGs in routine clinical practice, the following formula was then developed to calculate the risk score for each patient: Risk score= (−0.548 × expression IKBKB) + (1.3 × expression ST13) + (−0.744 × expression TSC2) + (−0.780 × expression MAP2K7). We found that the coefficient of IKBKB, TSC2, and MAP2K7 is negative, indicating that the expressions of IKBKB, TSC2, and MAP2K7 were related to better OS of laryngeal cancer patients. While the expression of ST13 was related to worse OS. Using the median expression value of risk score, the laryngeal cancer patients were divided into high-risk and low-risk groups.

The performance of the risk score was identified to predict the clinical outcome of laryngeal cancer patients, K–M plots were drawn to analyze the different survival rates between the high-risk and low-risk groups. The result of the K–M analyses indicated that patients in the high-risk group have a worse survival rate than those in the low-risk group (P = .005, Fig. 3C). Figures 3D–G showed the risk distribution of patients in the dataset, the number of patients in high- and low-risk groups, the survival time of patients in the dataset, and the heatmap of the 4 ARGs expression profiles in the dataset. Furthermore, univariate analyses and multivariate analyses were performed to assess the correlation between clinicopathological features and OS. After that, Figure 4A, B and Table 3 displayed that gender (HR = 0.401, 95% CI = 0.174–0.922, P = .031), N (HR = 1.631, 95% CI = 1.242–2.141, P < .001), and risk score (HR = 1.731, 95% CI = 1.312–2.280, P < .001) are very likely associated with patients OS.

Figure 4

Forest plots of clinical characteristics, receiver operating characteristic (ROC) curve of clinical characteristics, and box plots of the correlation between 4 key ARGs and clinical characteristics. (A) Univariate analyses of clinical characteristics. (B) Multivariate analyses of clinical characteristics. (C-G) Box plots with significant difference were displayed.

Table 3

Univariate and multivariate analyses of OS in laryngeal cancer patients of TCGA.

Furthermore, the correlations between clinicopathological features and ARGs expression were analyzed (Table 4). The analyses based on independent sample t tests indicated that the expression of IKBKB is higher in the elder patients (P = .010, Fig. 4C) and higher N (P = .054, Fig. 4D); the expression of MAP2K7 is higher in the higher stage group (P = .002, Fig. 4E) and higher N group (P = .005, Fig. 4F). Besides, the risk score values were higher in the higher stage group (P = .016, Fig. 4G).

Table 4

t tests of the correlation between 4 key ARGs and clinical characteristics in laryngeal cancer patients of TCGA.

Development of a nomogram for predicting prognosis of laryngeal cancer patients

We constructed a nomogram based on OS, which integrated the 4 ARGs (Fig. 5A). The nomogram can be interpreted by the total points assigned to each variable, which is indicated at the top of the scale. The total points can be converted to predict 1-year, 3-year, and 5-year OS for patients. The predictive accuracy of the nomogram is shown in Figure 5B. The AUC at 1-year prediction was 0.656, the AUC at 3-year prediction was 0.779, and the AUC at 3-year prediction was 0.796. Calibration curves for the nomogram revealed no deviations from the reference line and no need for recalibration (Fig. 5C (1-year), 5D (3-year), 5E (5-year)). Figure 5 F–H displayed the expression difference between high-risk and low-risk groups and the visualization of the correlation between risk score and survival time.

Figure 5

Construction of the nomogram, evaluation of autophagy-related genes (ARGs) prognostic model, and visualizations of risk score in the high-risk group and the low-risk group. (A) A nomogram of the prognostic model using 4 key ARGs. (B) Receiver operating characteristic (ROC) curve of 1-year, 3-year, and 5-year survival. (C-E) Calibration curves of 1-year, 3-year, and 5-year survival. (F-H) Visualizations of risk score in the high-risk group and the low-risk group.

Survival analyses of the autophagy-related model

Based on the TCGA database, we conducted the survival analyses of the autophagy-related model (IKBKB, ST13, TSC2, and MAP2K7). Then, Figure 6 showed that the overexpressing of IKBKB (P = .047) and MAP2K7 (P = .003) may indicate a better OS. However, there were no statistical differences in the survival analyses of ST13 (P = .268) and TSC2 (P = .094).

Figure 6

Survival analyses of 4 key autophagy-related genes (ARGs) based on The Cancer Genome Atlas (TCGA) database. (A)IKBKB. (B)ST13. (C)TSC2. (D)MAP2K7.

Discussion

Laryngeal cancer is a major malignancy of head and neck. The blank area in the knowledge of molecular targeted therapy and molecular biomarkers for laryngeal cancer prognosis monitoring still eagerly needed us to promote a better understanding of the molecular mechanisms.[ Exploration of the autophagy mechanism would open a new positive outlook for laryngeal cancer.[ However, most researches focused on autophagy only analyzed a signal gene or an autophagic inhibitor. Garcia-Mayea et al performed an immunohistochemistry study to identify the putative relevance of autophagy in laryngeal cancer.[ The expressions of the proteins (LC3, ATG5, p62/SQSTM1, and PTOV1) involved in the process of autophagy were analyzed. Finally, they suggested that autophagy inhibition with hydroxychloroquine could be a promising therapy for laryngeal cancer patients. Ji et al also indicated the protective role of autophagy in laryngeal cancer cells and provided a new attractive therapeutic strategy for laryngeal cancer with autophagic inhibitors.[ Also, Lin et al and Cao et al both proved the protective role of autophagy in laryngeal cancer.[

Given understanding advances in high-throughput sequencing, TCGA as the large-scale database provided effective measures for screening prognostic ARGs.[ In our study, we extracted the expression profiles of ARGs from TCGA and then aimed to select molecular biomarkers related to the prognosis of laryngeal cancer. Subsequently, we screened 38 differentially expressed ARGs between laryngeal cancer and non-tumor tissues. Considering these ARGs may be involved in the initiation of laryngeal cancer, GO and KEGG analyses of these genes were performed. Interestingly, the top 3 significant KEGG pathways (EGFR tyrosine kinase inhibitor resistance, apoptosis, and ErbB signaling pathway) of these enriched genes were decreased. Plenty of researches had proved that EGFR tyrosine kinase inhibitor resistance, apoptosis, and ErbB signaling pathway played key roles in various cancers.[ Thus, we indicated that autophagy may likely be considered as a tumor suppressor in the process of laryngeal cancer tumorigenesis. It is acknowledged that autophagy always inhibited tumorigenesis by activation of gene mutations or inactivation of ARGs.[

To obtain the autophagy-related molecular biomarkers for laryngeal cancer, we screened all ARGs and then identified 4 key prognostic ARGs, all of which may be the potential therapeutic targets. The result of the univariate survival analysis revealed that 4 ARGs were likely associated with OS in the TCGA database. Further multivariate survival analysis helped us to determine 4 key prognostic ARGs (IKBKB, ST13, TSC2, and MAP2K7) to develop the risk score, which could be an independent prognostic indicator for laryngeal cancer patients. Furthermore, we further analyze prognostic ARGs and clinical characteristics, which would provide an accurate estimation of OS in laryngeal cancer. Thus, univariate analyses and multivariate analyses based on clinicopathological features showed that gender, N, and risk score are very likely associated with patients OS. The analyses based on independent sample t tests showed that the expression of IKBKB is higher in the elder patients and higher N; the expression of MAP2K7 is higher in the higher stage group and higher N group. In addition, risk score values were higher in the higher stage group. To help doctors and patients to perform an individualized survival prediction, a nomogram was constructed to predict the 1-, 3-, and 5-year OS by using this 4 key ARGs and clinicopathological risk variables. These findings could also provide an effective multi-dimensional biomarker strategy, which would be effective in monitoring autophagy and predicting the prognosis in laryngeal cancer patients. Last, we conducted the survival analyses of significant prognostic ARGs (IKBKB, ST13, TSC2, and MAP2K7). Based on the TGGA database, we found that the overexpressing of IKBKB and MAP2K7 may indicate a better OS. However, there were no statistical differences in the survival analyses of ST13 and TSC2. Therefore, the prognostic signature could be converted into a clinical application.

IKBKB, inhibitor of nuclear factor-kappa b kinase subunit beta, is a protein-coding gene that played an essential role in the NF-kappa-B signaling pathway. Krazinski BE et al indicated that IKBKB protein could be associated with clinical relevance in clear cell renal cell cancer, acting as a marker of poor prognosis.[ Li et al proved that tumor necrosis factor significantly increased phosphorylation of IKBKB.[ MAP2K7, mitogen-activated protein kinase kinase 7, a protein-coding gene that acts as an essential component of the MAP kinase signal transduction pathway. Shen et al uncovered that MAP2K7 was repressed by KLF4 in T-cell acute lymphoblastic leukemia cells.[ Hong et al and Min et al proved that MAP2K7 may be a potential target for tumor therapy.[ ST13, st13 hsp70 interacting protein, is a protein-coding gene that suggested that it may be a candidate tumor suppressor gene. Yu et al, Yang et al, and Wang et al all concluded that ST13 may be a novel therapeutic target for colorectal cancer.[ TSC2, TSC Complex Subunit 2, is a protein-coding gene that is related to the pathways of Vesicle-mediated transport and mTOR signaling.[ Schrader et al and Li et al proved that TSC2 could be a tumor suppressor gene.[ Hence, combined with our results, we concluded that IKBKB, ST13, TSC2, and MAP2K7 may play the tumor suppressor roles, and IKBKB and MAP2K7 may be the most vital.

However, there are indeed several potential limitations in our study. Firstly, owing to the lack of sufficient cases and reliable laryngeal cancer cells (the standard laryngeal cancer cell - Hep-2 cell was challenged), we failed to validate the expression of IKBKB, ST13, TSC2, and MAP2K7. Secondly, our study was based on the data from the TCGA database. Thirdly, the information on several other important clinical outcomes, such as various treatments, and the number of lymph nodes, was unavailable now.

Based on these comprehensive analyses with ARGs expression profiles and clinical outcomes, 4 prognostic ARGs (IKBKB, ST13, TSC2, and MAP2K7) were finally identified. A novel autophagy-related model and a predictive nomogram were conducted to robustly estimate laryngeal cancer patients survival. However, further experiments are still needed to test our conclusion. Supplementary Table 1.

Author contributions

Conceptualization: Meng-Si Luo, Guan-Jiang Huang, Hong-Bing Liu.

Data curation: Meng-Si Luo, Guan-Jiang Huang.

Formal analysis: Meng-Si Luo, Guan-Jiang Huang, Hong-Bing Liu.

Investigation: Meng-Si Luo, Guan-Jiang Huang, Hong-Bing Liu.

Methodology: Meng-Si Luo, Guan-Jiang Huang, Hong-Bing Liu.

Project administration: Hong-Bing Liu.

Resources: Meng-Si Luo, Guan-Jiang Huang, Hong-Bing Liu.

Software: Meng-Si Luo, Guan-Jiang Huang.

Supervision: Meng-Si Luo, Hong-Bing Liu.

Validation: Meng-Si Luo.

Visualization: Guan-Jiang Huang.

Writing – original draft: Meng-Si Luo, Guan-Jiang Huang.

Writing - review & editing: Meng-Si Luo, Guan-Jiang Huang, Hong-Bing Liu.