Differentially expressed ARGs between prostate cancer and adjacent non-tumor tissues
A total of 485 primary PCa patients with RNA-seq data and clinical follow-up information were involved in the present study. Among 234 autophagy-related genes, there were 13 differentially expressed ARGs, including 5 up-regulated (ATG9B, BIRC5, CAMKK2, CDKN2A, and NKX2-3) and 8 down-regulated ARGs (DNAJB1, FAM215A, HSPB8, ITGB4, ITPR1, NRG1, NRG2, and TP63), with thresholds of |log2 fold change (FC)| > 2 (Fig. 1a). Then, the volcano plot and box plot were visualized to show the expression pattern of the differentially expressed ARGs between PCa and non-tumor tissues (Fig. 1b, c).
GO enrichment analysis of differentially expressed autophagy-related genes
GO enrichment analysis was performed according to the differentially expressed ARGs. According to the results of DAVID, we found that the top enriched GO terms for biological processes were autophagy, process utilizing autophagic mechanism, and odontogenesis of dentin-containing tooth. The heatmap of the relationship between ARGs and GO enrichment analysis was also displayed (Fig. 2a, b).
Identification of prognosis-related ARGs and construction of prognosis prediction model
A total of 14 ARGs were significantly associated with OS in the univariate Cox regression analysis. Furthermore, in the multivariate Cox regression analysis, five genes including FAM215A, FDD, MYC, RHEB, and ATG16L1 were identified to construct the OS prediction model. OS-related prediction model = (17.20896* expression value of FAM215A) + (4.319028* expression value of FADD) + (0.674838* expression value of MYC) + (1.869633* expression value of RHEB) + (2.071004* expression value of ATG16L1).
We divided the 485 prostate cancer cases into high- and low-risk groups according to the median values of the OS-related prediction model. Kaplan–Meier survival curves showed that low-risk group had a lower mortality rate than high-risk group (HR = 6.391, 95% CI = 1.581–25.840, P < 0.001) (Fig. 3a). The ROC curves of OS-related predictive signatures were demonstrated in Fig. 3b, with AUC of 0.84. Figure 3c, d showed the OS-related prediction model distribution of patients in the TCGA dataset.
According to the median value of the five genes, the high expression level of FAM215A (HR = 4.347, 95% CI = 1.175–16.290, P = 0.041), FADD (HR = 7.009, 95% CI = 1.892–25.960, P = 0.031), and MYC (HR = 7.153, 95% CI = 1.932–26.470, P = 0.029) were significantly associated with worse OS in Kaplan–Meier curves (Fig. 4). However, this association did not hold true of gene ATG16L1(HR = 2.426, 95% CI = 0.653–9.017, P = 0.194) and RHEB (HR = 1.236, 95% CI = 0.335–4.566, P = 0.744) in Kaplan–Meier curves (Additional file 1: Fig. S1).
Among 234 autophagy-related genes, a total of 53 ARGs were significantly associated with DFS in the univariate Cox regression analysis. In the multivariate Cox regression analysis, a total of 22 genes were significantly associated with DFS in PCa (Fig. 5a). DFS-related prediction model = (0.97225* expression value of ULK2) + (− 1.74297* expression value of NLRC4) + (− 1.11799* expression value of MAPK1) + (− 1.12182* expression value of ATG4D) + (− 0.73348* expression value of MAPK3) + (1.40252* expression value of ATG2A) + (− 0.49364* expression value of ATG9B) + (− 1.09886* expression value of FOXO1) + (− 0.68955* expression value of PTEN) + (1.80095* expression value of HDAC6) + (− 0.99993* expression value of PRKN) + (0.35846* expression value of HSPB8) + (− 0.51552* expression value of P4HB) + (1.56551* expression value of MAP2K7) + (− 0.96348* expression value of MTOR) + (1.65516* expression value of RHEB) + (0.73934* expression value of TSC1) + (0.27799* expression value of BIRC5) + (1.43484* expression value of RGS19) + (− 0.63037* expression value of RAB24) + (− 0.28580* expression value of PTK6) + (− 1.05312* expression value of NRG2).
We divided the PCa cases into high- and low-risk groups according to the median values of the DFS-related prediction model. Kaplan–Meier survival curves showed that high-risk group had a lower disease-free rate than low-risk group (HR = 7.407, 95% CI = 4.850–11.320, P < 0.001). The ROC curves of OS-related predictive signatures were demonstrated in Fig. 5b, with AUC of 0.85. Figure 5c, d showed the DFS-related prediction model distribution of patients in the TCGA dataset.
Among the 22 genes in DFS-related prediction model, high expression of ATG2A (HR = 2.266, 95% CI = 1.492–3.442, P < 0.001), ATG4D (HR = 1.665, 95% CI = 1.096–2.530, P = 0.017), ATG9B (HR = 1.803, 95% CI = 1.187–2.738, P = 0.007), BIRC5 (HR = 2.013, 95% CI = 1.384–3.195, P < 0.001), MAPK3 (HR = 2.148, 95% CI = 1.414–3.263, P < 0.001), NLRC4 (HR = 2.053, 95% CI = 1.352–3.119, P = 0.001), RAB24 (HR = 2.811, 95% CI = 1.851–4.270, P < 0.001), RGS19 (HR = 2.019, 95% CI = 1.329–3.068, P = 0.001), RHEB (HR = 2.137, 95% CI = 1.407–3.245, P < 0.001), ULK2 (HR = 1.579, 95% CI = 1.039–2.399, P = 0.033), and TSC1 (HR = 1.622, 95% CI = 1.067–2.464, P = 0.024) genes were associated with worse prognosis in PCa in Kaplan–Meier curves according to the median values of gene expression (Fig. 6). In addition, high expression of FOXO1 (HR = 2.087, 95% CI = 1.373–3.172, P < 0.001), HSPB8 (HR = 1.673, 95% CI = 1.101–2.541, P = 0.017), MTOR (HR = 1.897, 95% CI = 1.247–2.885, P = 0.002), NRG2 (HR = 1.944, 95% CI = 1.280–2.955, P = 0.002) and PRKN (HR = 2.308, 95% CI = 1.518–3.508, P < 0.001) genes were associated with better prognosis in Kaplan–Meier curves according to the median values of gene expression (Fig. 7). No differences were found between the expression level of HDAC6 (HR = 1.392, 95% CI = 0.913–2.123, P = 0.116), MAP2K7 (HR = 1.379, 95% CI = 0.908–2.094, P = 0.133), MAPK1 (HR = 1.426, 95% CI = 0.939–2.167, P = 0.095), P4HB (HR = 1.501, 95% CI = 0.988–2.280, P = 0.058), PTK6 (HR = 1.338, 95% CI = 0.881–2.032, P = 0.174), and PTEN (HR = 1.324, 95% CI = 0.872–2.010, P = 0.191) and disease-free survival (Additional file 2: Fig. S2).
The relationships between clinicopathological parameters and prognosis-related ARGs and prognosis-related prediction model
The OS-related prediction model values were higher in T3-4 than in T1-2 (P = 0.008), and higher in Gleason score > 7 than ≤ 7 (P = 0.015). No difference of OS-related prediction model values was observed between age > 65 than age ≤ 65 (P = 0.164), or N0 stage and N1 stage (P = 0.088) (Fig. 8). The DFS-related prediction model values were higher in T3-4 than in T1-2 (P < 0.001), higher in N1 than in N0 (P = 0.001), and higher in Gleason score > 7 than ≤ 7 (P < 0.001). No difference of DFS-related prediction model values was observed between age > 65 than age ≤ 65 (P = 0.208) (Fig. 9).
Among 485 primary PCa patients in the present study, only two of them had distant metastasis. Therefore, the relationship between the M status and prediction model has not been analyzed.