Identification of an autophagy-related gene expression signature for colorectal cancer

Abstract

Autophagy, a major lysosomal degradation process, has been implicated in the pathogenesis of colorectal cancer (CRC). However, it still remains unclear what role all the autophagy-related genes (ARGs) play in CRC prognosis. A univariate Cox analysis was conducted to select ARGs related with overall survival (OS) for CRC cohort of the TCGA database. A multivariate Cox analysis was used to develop an ARGs expression signature, whose risk score stratified CRC patients into two groups. The OS time between the above two groups was estimated using the Kaplan–Meier curve and log-rank test. In order to evaluate the predictive performance of this expression signature and other clinicopathological parameters for CRC prognosis, the area under the curve (AUC) was calculated. Nine out of 20 OS-related ARGs were finally incorporated into this ARGs expression signature, whose risk score stratified CRC patients into either high-risk or low-risk group. The OS time of patients in low-risk group was higher than that in high-risk group and the difference was statistically significant (P < 0.001). The AUC for this risk score was higher than other clinicopathological parameters, including TNM stage. This study presents an ARGs expression signature, which outperforms other clinicopathological parameters in predicting CRC prognosis.

KEYWORDS:

• Autophagy
• colorectal cancer
• prognostic signature
• survival analysis
• TCGA

Introduction

Autophagy is a major catabolic mechanism by which denatured, damaged or aging organelle and protein are delivered to lysosome for self-degradation in eukaryotic cells (Bagherniya et al. 2018; Packer 2020). By eliminating these abnormal macromolecular components, autophagy contributes to the maintenance of intracellular homeostasis (Fang et al. 2020). Autophagy abnormality participates in the development of many diseases, including immune dysfunction, neurodegenerative disease, chronic disease and tumor (Ait-Aissa et al. 2020; Li et al. 2020; Zhao and Zhang 2019). However, it still remains controversial what role autophagy plays in tumor (Levy et al. 2017; Xu et al. 2020). In premalignant stage, stimulation of autophagy could prevent cancer development. However, increased autophagic flux enables tumor growth in advanced cancer. Thus, it is valuable to further explore the role of autophagy in tumor, which could help provide the important lessons for cancer therapies by deliberately manipulating autophagy.

Colorectal cancer (CRC) is currently one of the most frequent malignant tumors (Abu et al. 2020; Mauri et al. 2019; Zhou et al. 2020). With the increasing incidence of CRC in the young, there will be an estimated 2.5 million newly diagnosed CRC cases in 2035. Although the current treatment options, including surgery, chemotherapy and molecular targeting agents have been used to extend survival for CRC patients, poor prognosis still remains a serious problem (Al Bandar and Kim 2017; Dekker et al. 2019). Although the TNM classification was frequently used to predict CRC prognosis (Amin et al. 2017), it could not sufficiently assess individual patient outcome. Moreover, patients within the same TNM stage could have a significant difference in clinical prognostication (Huang et al. 2019). Therefore, it is essential to identify novel potential biomarkers and gene expression signature for CRC prognosis.

Recently, accumulating evidence has reported that autophagy is implicated in CRC development and progression, which could provide a new insight for CRC therapeutic strategies by selectively targeting autophagy (Mokarram et al. 2017). For example, Wei et al. found that FAT4 could suppress epithelial–mesenchymal transition (EMT) of CRC by stimulating autophagy in part through the P13K-AKT signaling pathway (Wei et al. 2019). Yu et al. revealed that gut microbiota could promote CRC resistance to chemotherapy by activating the autophagy pathway (Yu et al. 2017). However, previous researches only highlighted the role of one or a few autophagy-related genes (ARGs) in CRC. Thus, the purpose of this research focuses on the development of a novel autophagy-related expression signature for CRC prognosis, based on a comprehensive analysis of the ARGs.

Materials and methods

Data acquisition

Two hundred and thirty-two human genes involved in autophagy were collected from the Human Autophagy Database (HADb, http://www.autophagy.lu/clustering/index.html). Both expression profiles of ARGs in RNA-sequencing (RNA-seq) and clinical data of CRC cohort were obtained and extracted from TCGA website (https://portal.gdc.cancer.gov/).

Enrichment analysis of differentially expressed ARGs

Differentially expressed ARGs between CRC and non-tumor samples were determined using limma package in R software. The adjust P-value < 0.05 and |log2 fold change (FC)| > 1 were chosen as the threshold criteria. Then, the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of these differentially expressed genes were performed to find the biological significance and essential pathways by clusterprofiler R package. To provide visual information, the enriched results for GO and KEGG analysis were presented as the bubble plot and circle map, respectively.

Development of an autophagy-related gene signature

Log2 (count+1) transformation was used to normalize ARGs expression profiles of TCGA-CRC cohort. Univariate Cox analysis was first conducted to select the ARGs remarkedly related with overall survival (OS) of TCGA-CRC cohort. Subsequently, multivariate Cox regression analysis was performed using the above OS-related ARGs. Some genes, which might not serve as an independent predictor in CRC prognosis, were removed. Finally, an autophagy-related gene signature incorporating several prognostic genes was developed based on ARGs expression values multiplied by Cox regression coefficients. CRC patients were then split into two groups with a cutoff point of the median risk score. The OS time between the above two groups was estimated using the Kaplan–Meier curve and log-rank test. To explore whether the risk score could independently predict OS of TCGA-CRC cohort, we performed the multivariate Cox analysis by survival R package. Moreover, receiver operating characteristic (ROC) curve with multiple clinicopathological variables were plotted and area under the curve (AUC) was measured to evaluate the predictive performance by survivalROC R package. Finally, this ARGs expression signature was transformed into a visual nomogram by rms R package.

Statistical analysis

R software was used to perform all the statistics. Statistical difference was found with P-value no more than 0.05. Univariate and multivariate Cox analysis were performed to screen the OS-related ARGs and develop the ARGs expression signature, respectively.

Results

Identification of differentially expressed ARGs

TCGA-CRC cohort consisted 568 tumor and 44 non-tumor samples, whose RNA-seq and clinical data were downloaded. A total of 541 CRC patients with available gene expression value and OS time were finally included in this research. Expression values of 232 ARGs were extracted. Thirty-seven differentially expressed ARGs were obtained with the threshold criteria of FDR < 0.05 and |log2 FC| > 1 (Figure 1(A)). To visualize the expression profiles of these genes between CRC and non-tumor samples, a box plot was displayed, in which 19 up-regulated genes (FAS, CAPN10, SESN2, PINK1, GABARAP, TP53INP2, ERO1A, CAPN2, TNFSF10, VEGFA, BID, BIRC5, HSP90AB1, MYC, CDKN1A, CD46, EIF4EBP1, ATIC, BCL2L1) and 18 down-regulated genes (GRID2, NRG3, MAP1LC3C, TMEM74, NRG2, IFNG, CCR2, NRG1, PRKN, TP73, ATG9B, CDKN2A, SPHK1, HSPB8, NKX2−3, FKBP1B, ITPR1, BCL2) were identified (Figure 1(B)).

Figure 1. Gene enrichment analysis of differentially ARGs. (A) Heat map shows differentially expressed ARGs between CRC and non-tumor samples. (B) Box plot identifies 19 up-regulated genes and 18 down-regulated genes. (C) GO analysis shows the biological process, cellular component, and molecular function of differentially expressed ARGs. (D) KEGG analysis shows the signaling pathway involved in differentially expressed ARGs.

Enrichment analysis of differentially expressed ARGs

To further understand the biological significance of 37 differentially expressed ARGs, we performed GO and KEGG functional enrichment analysis. As was shown in Figure 1(C), the top GO enrichment terms for molecular function were: ubiquitin protein or ubiquitin-like protein ligase binding, kinase or protein kinase regulator activity; for cellular components were: autophagosome, vacuolar membrane, and autophagosome membrane; and for biological process were: autophagy, process utilizing autophagic mechanism and response to oxygen levels. Besides, the KEGG results revealed that these genes were mainly enriched in p53 signaling pathway, platinum drug resistance, apoptosis and the like (Figure 1(D)).

Identification of prognostic ARGs

The association between ARGs expression value of 541 TCGA-CRC patients and OS was investigated, resulting in 20 OS-related ARGs. As was shown in Figure 2(A), most of these 20 genes were identified as risk factors. Considering the potential clinical significance, we further investigated the genetic mutations of these ARGs via cBioportal database. The results indicated that among 592 CRC patients, genetic mutation was found in 483 patients and the mutation rate was as high as 82%. The top five genes with genetic mutation were WIPI2, ATIC, RAB7A, NRG1 and MAPK9, with mutation rate of 35, 16, 15, 14 and 12%, respectively (Figure 2(B)).

Figure 2. Genetic alteration of OS-related ARGs. (A) Forest plot shows 20 OS-related ARGs of TCGA-CRC cohort by univariate Cox regression analysis. (B) Top 5 ARGs with genetic mutation were WIPI2, ATIC, RAB7A, NRG1 and MAPK9, with mutation rate of 35, 16, 15, 14 and 12%, respectively.

Based on the above prognostic ARGs in univariate analysis, multivariate analysis was conducted to develop the ARGs expression signature, whose risk score was [DAPK1 expression * (0.435964094)] + [ULK3 expression * (0.395252108)] + [ATIC expression * (0.630860629)] + [NFKB1 expression * (−0.584152048)] + [RAB7A expression * (0.982036429)] + [WIPI1 expression * (−0.261064362)] + [EDEM1 expression * (−0.626851909)] + [CX3CL1 expression * (0.175440139)] + [GRID1 expression * (1.48676544)]. By calculating the risk value of each patient, we divided CRC patients into two groups based on the median risk score. As was shown in Figure 3(A), the patients in low-risk group had longer OS time than that in high-risk group. Besides, the heatmap of prognostic ARGs expression profiles in multivariate analysis, distribution of risk score of CRC patients and relationship between survival time and risk score were presented in Figure 3(B,C). The AUC for risk score was 0.716, which was higher than tumor stage (AUC_stage = 0.709), TNM stage (AUC_T= 0.673, AUC_N= 0.656, AUC_M= 0.647), age (AUC_age= 0.646) and gender (AUC_gender= 0.433) (Figure 3(D)). Then, the value of the risk score in prognosis of CRC was evaluated. Univariate analysis revealed that the risk score was remarkedly related with OS (Figure 4(A)). Multivariate analysis showed that the risk score could independently predict CRC prognosis (Figure 4(B)).

Figure 3. Characteristics of this autophagy-related gene signature. (A) Kaplan–Meier analysis reveals that patients in low-risk group have longer OS time than that in high-risk group. (B) Heatmap of the expression profile of prognosis-related ARGs. (C) Distribution of risk score of CRC patients and relationship between survival time and risk score. (D) The AUC for risk score is higher than other clinicopathological parameters, including TNM stage.

Figure 4. Predictive performance of this autophagy-related gene signature. A forest plot in univariate (A) and multivariate(B) Cox regression analysis shows that this expression signature could serve as an independent predictor for CRC prognosis.

Clinical utility of prognostic signature and construction of the nomogram

Finally, the association between risk score and various clinicopathological parameters was analyzed. Results revealed that the risk score was positively related with TNM stage (P_T= 0.025, P_N < 0.001, P_M= 0.009) and tumor stage (P < 0.001) (Figure 5). To help clinicians predict prognosis of CRC patients more conveniently, the ARGs expression signature was finally presented as a visual nomogram (Figure 6).

Figure 5. Clinicopathological significance of this autophagy-related gene signature in CRC. This risk score is positively related with tumor stage, tumor size, lymph node involvement and distant metastasis.

Figure 6. Development of a nomogram to predict the risk of OS in CRC.

Discussion

CRC is a major lethal malignant tumor that urgently requires effective prognostic biomarkers (Shi et al. 2018). Accumulating evidence has identified the role of autophagy in the pathogenesis of CRC, which opens a new perspective for its diagnosis and prognosis (Wei et al. 2019; Yu et al. 2017). However, most of the previous studies only focused on the role of a single gene from the perspective of autophagy. A comprehensive analysis of relationship between all the ARGs expression profiles and CRC prognosis was not performed.

With the advances in sequencing technology, many databases, including TCGA and GEO, were used to screen significant autophagy-related gene signatures (An et al. 2018; Mo et al. 2019; Wang et al. 2019). To search molecular biomarkers for CRC prognosis, the ARGs expression profiles of TCGA-CRC cohort were deeply mined in our research. First, by comparing expression value between CRC and non-tumor samples, 37 differentially expressed ARGs were identified, of which GO and KEGG analysis were performed. Interestingly, the enrichment results revealed that ubiquitin protein ligase binding and p53 signaling pathway were the most significant molecular function and KEGG pathway. On the basis of these results, we hypothesized that these differentially ARGs may function as ubiquitin protein ligase, which could bind to p53 tumor suppressor gene, resulting in its ubiquitylation and subsequent degradation (Khoo et al. 2014). Dysregulation of p53 was involved in the pathogenesis of CRC (Li et al. 2015). This may provide a new insight into targeting p53-binding ARGs for CRC therapy.

Twenty ARGs were found associated with OS of TCGA-CRC cohort in univariate survival analysis. Then, nine key prognosis-related ARGs identified in multivariate analysis were used to develop an autophagy-related gene signature, which could serve as an independent predictor in CRC prognosis. Besides, this gene signature was useful for stratifying outcome of CRC patients. Accumulating evidence revealed that the existing TNM stage could not accurately predict the prognosis between stage II c and III a in CRC, partly due to overestimating the weight of N stage and underestimating the role of T stage (Li et al. 2014; 2016). Our research revealed that this ARGs expression signature outperformed other clinicopathological parameters including TNM stage, in predicting survival of CRC patients, which could provide a new strategy for CRC prognosis.

Most of these nine genes have been found related with the prognosis of CRC or other malignant tumors. DAPK1 was reported to be a metastasis suppressor, which was associated with its role of inhibiting integrin-mediated cell adhesion and migration. Decreased DAPK1 protein was found related with the metastasis and worse prognosis of multiple malignant tumors including CRC (Steinmann et al. 2019). ULK3, as a serine and threonine kinase, was involved in the phosphorylation of GL1 proteins to promote their transcriptional activity, resulting in positive regulation of Hedgehog pathway (Maloverjan et al. 2010). Abnormal activation of Hedgehog pathway was strongly implicated in the pathogenesis of various cancers, including CRC (Chowdhury et al. 2013). ATIC, as a key enzyme involved in the de novo purine synthesis, was strongly related with the tumor proliferation. Li et al reported aberrant ATIC upregulation would result in proliferation, migration and colony formation of hepatocellular cells by affecting the AMPK-mTOR-S6 K1 pathway, a pivotal axis protecting against various pathological alterations (Li et al. 2017; 2020). NFKB1 worked as a transcription factor, which was strongly implicated in the relationship between chronic inflammation and CRC (Patel et al. 2018). Cooks et al. demonstrated NF-κB activation induced p53 suppressor gene mutation, which would aggravate inflammatory response and participate in inflammation-related CRC carcinogenesis (Cooks et al. 2013). RAB7A, as a member of RAS oncogene family, encoded Rab-7a protein, which participated in the multiple biological processes, including apoptosis, cell cycle arrest and EMT. Xie et al. revealed that Rab7a overexpression would promote the ability of proliferation, growth and migration in breast cancer (Xie et al. 2019). CX3CL1 was a multifunctional chemokine, which could promote cancer survival and recruit anti-tumor leukocytes. In CRC, Conroy et al. reported higher CX3CL1 was positively related with higher tumor infiltrations of lymphocytes, resulting in inhibition of metastasis and better prognosis (Conroy and Lysaght 2020). Su et al. reported that no relationship was found between the polymorphism of GRID1 and CRC in the Chinese Han population (Su et al. 2015). However, in our research, GRID1 overexpression was found significantly related with advanced TNM stage and tumor stage. Besides, GRID1 could serve as an independent prognostic factor. Thus, future researches were required to explore the role of GRID1 in CRC.

Several shortcomings should be acknowledged in our study. The most obvious is that some inherent bias may be introduced by the retrospective nature. Second, the potential mechanisms between nine prognosis-related genes and CRC are not elucidated. Future studies are required to explore the underlying mechanisms. Finally, this autophagy-related gene signature is not validated in other CRC cohorts. Further researches are needed to ensure the robustness of our ARGs expression signature.

Conclusion

In summary, based on TCGA-CRC cohort, we develop an expression signature based on nine ARGs, which is an independent prognostic indicator for CRC prognosis. Additionally, this ARGs signature could stratify CRC patient outcome and outperform other clinicopathological parameters including TNM stage. These ARGs signature may help accelerate personalized treatment. However, future prospective researches are required to assess its clinical utility.

Acknowledgements

All authors thank the family for participating in the study. No sources of funding were required for this work.

Authors’ contributions

L.X. performed the statistical procedure and drafted the paper, Y.Q. helped edit the paper, and Q.Z. provided language help.

Disclosure of interest

The authors report no conflict of interest.

Data Availability Statement

Data sharing is not applicable to this article as no new data were created or analyzed in this study.