Distinct autophagy-associated genes in EAC and normal tissues
A total of 9 normal and 78 EAC tissues with gene expression profiles and clinical data were obtained from TCGA. There were 30 significantly different autophagy-related genes (SD-ARGs) between the normal and tumor groups. Among these genes, 6 genes were down-regulated and 24 were up-regulated in the tumor group compared with normal group (Table 1). The heatmap, volcano plot and bar plot were shown in Figure 1a-c.
Table 1 about here
Figure 1 about here
Enrichment analysis
GO (Gene Ontology) analysis included the biological process (BP), cellular component (CC) and molecular function (MF) categories. In the BP category, the SD-ARGs were obviously enriched in the intrinsic apoptotic signaling pathway, as well as the regulation of apoptotic signaling pathway. In the CC, the SD-ARGs were mainly enriched in the mitochondrial outer membrane process. In addition, the CC analysis also showed the SD-ARGs are involved in the autophagosome and autophagosome membrane. In the MF, the SD-ARGs were obviously enriched in ubiquitin protein ligase binding process (Figure 2a). In the KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis, the SD-ARGs were mainly enriched in the apoptosis, which was similar in the GO analysis (Figure 2b).
In the GO circle plot, the SD-ARGs were mainly enriched in the intrinsic apoptotic signaling pathway and the regulation of cytochrome c from mitochondrion. And in the KEGG circle plot, it was also clearly showed these genes were mainly involved in the apoptotic process. Many SD-ARGs have other pathways related to the autophagy, such as the apoptosis-multiple species shown in the KEGG circle plot (Figure 2c-d).
The GO and KEGG heatmaps revealed the SD-ARGs were enriched in the apoptotic process (Figure 2e-f).
Figure 2 about here
Prognosis-related ARGs
Univariate cox regression was used to investigate ARGs related with prognosis, 14 ARGs (TBK1, CAPN2, ATG5, GOPC, TP73, BECN1, RB1CC1, SIRT1, VAMP7, DDIT3, DAPK1, ATG12, CAPN1, ITGA3) were found to be significantly associated with overall survival (OS) (shown in Figure 3a). Next, multivariate cox regression was performed to select the proper ARGs from these 14 genes. After calculating, eight ARGs (ATG5, TP73, BECN1, SIRT1, VAMP7, DAPK1, ATG12, CAPN1) was selected (shown in Table 2). Among these eight genes, BECN1, DAPK1 and CAPN1 played protective roles (HR<1), while the other five genes were risk factors in survival (HR>1). According to the eight genes expressions and their coefficient [18], we then calculated the risk score (= , with Coef j indicating the coefficient and Xj representing the relative expression levels of each ARG standardized by z-score) of each patient and used the median risk score value as a cutoff point for classifying the EAC patients into a high-risk group (n = 39) and a low-risk group (n = 39), respectively. Patients in the high-risk group obviously had a shorter overall survival time than patients in the low-risk group (median time = 1.10 years vs. 1.752 years, p < 0.001, Fig. 3b).
Figure 3 about here
Table 2 about here
Prognostic hazard curves
All the tumor patients were divided into high-risk group or low-risk group. As the risk score increased, the patients’ risk increased, and the survival time decreased (Figure 4a-b). The risk heatmap clearly showed CAPN1 was down-regulated in high-risk group, implying a tumor-suppressor role. However, VAMP7 was up-regulated in high -risk group and this implied it was a tumor-promoting role.
Figure 4 about here
Independent risk factors of OS
We combined the risk score with clinical data in EAC patients. Then we performed univariate cox and multivariate cox regression analyses to investigate the independent risk factors for OS. The univariate cox analysis showed that tumor stage, M (metastasis), N (lymph nodes) and risk score were correlated with OS (P<0.001, P=0.002, P=0.047, P<0.001, respectively) (Figure 5a). Multivariate cox regression showed gender (HR=0.225, P=0.032), stage (HR=5.841, P=0.008) and risk score (HR=1.131, P<0.001) were independent risk factors for survival (Figure 5b).
Figure 5 about here
Model development for predicting the survival
In order to provide an approach to predict the survival, we constructed the ROC curve using the independent risk factors associated with OS (gender, stage and risk score). In addition, we assess the feasibility using the area under curve (AUC) values. The risk score curve showed better ability to predict the survival risk (AUC=0.892) than other indicators (Figure 6).
Figure 6 about here
Clinical correlation analysis
In order to explore the relationships between the prognostic eight ARGs and clinical features, we calculated the correlations using the t-test or Kruskal-Wallis test. We found BECN1, DAPK1, VAMP7 and risk score were significantly associated with survival status (all P values <0.05). Among these, BECN1 and DAPK1 expressions were higher in the alive patients, implying their protective roles (Figure 7a-d). We also found SIRT1 was significantly associated with gender, tumor stage and T (primary tumor) and its expression levels increased with time (Figure 7e-g). In addition, risk score was also associated with tumor stage (Figure 7h).
Figure 7 about here
Experimental validation
According to the screening and validation steps as described above, we selected the four most significant ARGs (BIRC5, CDKN2A, ITPR1,PRKN) from the 30 significantly different genes according to the P (<0.0001) and FDR values (<0.0001), and GAPDH was set as an internal reference. The experiment results showed the BIRC5 was overexpressed in the EAC tissues and the ITPR1 and PRKN were downregulated in the normal esophageal mucosal tissues. These results were shown in Figure 8a-c.
Figure 8 about here