3.1 ATG signature establishment
After filtration with MAD >0.5, 617 genes measured on all platforms were selected for this study. By 1000 times of random Cox univariate regressions, 58 ATGs were identified to be strongly relevant to DFS and considered as candidates for the signature (Figure 1). A LASSO Cox regression in stage I and II patients in the training cohort (Table 1) revealed 10 ATGs for the risk model (Figure 2), with coefficient of each ATG listed in Table 2. The risk model was formulated as follow: Autophagy-related risk score = 0.040346631× exp (mRNA expression level of CD163L1) + 0.040346631× exp (mRNA expression level of FAM13B) + 0.160103165× exp (mRNA expression level of HDAC6) + 0.05063732× exp (mRNA expression level of HPR) − 0.012205947× exp (mRNA expression level of NR2C2) − 0.027104325× exp (mRNA expression level of RAB12) + 0.055095935× exp (mRNA expression level of SIRT2) − 0.084226324× exp (mRNA expression level of TBC1D14) − 0.012613054× exp (mRNA expression level of TLK2) − 0.001301898× exp (mRNA expression level of TBC1D12).
Based on the time-dependent ROC curve analysis at 5 year-DFS in the training cohort, the optimal cutoff of ATG signature that divided the patients into high and low autophagy risk groups was defined as -0.009 (Figure 3). The risk scores of all patients are shown in Supplementary Table 1. Survival analysis showed the DFS rate was higher in the low autophagy risk group compared to the high autophagy risk group for patients with early stages (stage I and II) in the training cohort (GSE39582) (Figure 4A, 4B, 4C, HR, 2.76[1.56–4.82]; P=2.06×10-4). So was indicated for patients in all stages in the GSE39582 dataset (Figure 5A, 5B, 5C, HR, 1.7[1.25–2.31]; P=5.21×10-4).
3.2 Validation of the ATG signature
To assess the predictive capability of the risk model, C-index was first applied to various cohorts which turned out to be 0.74 (95%CI, 0.63 - 0.85) in the GSE39582 cohort, 0.70 (95%CI, 0.54 - 0.85) in the TCGA cohort and 0.70 (95%CI, 0.51 - 0.89) in the meta-validation cohort (Table 3), higher than those of Oncotype DX colon. To further verify whether the ATG signature predicted well in various populations, we employed the same formula to the independent validation cohort (TCGA) and the meta-validation cohort (GSE37892 and GSE14333). Patients were significantly stratified into different risk groups by the ATG signature considering DFS. For early stages, CRC patients in high autophagy risk group displayed poorer DFS in both the independent validation cohort (Figure 4D, 4E, 4F, HR, 2.29[1.15–4.55]; P=1.5×10-2) and the meta-validation cohort (Figure 4G, 4H, 4I, HR, 2.5[1.03–6.06]; P=3.63×10-2). So was indicated for patients with all stages in both the independent validation (Figure 5D, 5E, 5F, HR, 1.79[1.16–2.7]; P=5.3×10-3) and the meta-validation cohorts (Figure 5G, 5H, 5I, HR, 1.64[1.04–2.52]; P=3.16×10-2). Besides, univariate and multivariate analysis further proved ATG as an independent prognostic factor after adjusting for clinical parameters such as sex and tumor stage (Table 4).
3.3 Functional analysis of the ATG signature
GO analysis and GSEA were carried out to explore the biological function and signaling pathways of genes from the ATG signature. GO analysis revealed that the genes within the ATG signature were mostly involved in the regulation of autophagy and catabolic processes (Figure 6A and Supplementary Table 2). GSEA was performed between different risk groups to further investigate the pathways that were significantly affected. We found a significant enrichment in multiple immune/ inflammatory pathways in the high autophagy risk group, including the IL6/JAK/STAT3 signaling pathway, the IL2/STAT5 pathway, the IFN-αpathway, the IFN-γ pathway, the inflammatory response pathway (Figure 6B, P value<0.005). Some cell cycle/ metabolism associated pathways, including the G2-M, oxidative phosphorylation, E2F, MYC were also significantly enriched in the high autophagy risk group, in addition to a few classic pathways like mTORC1 and epithelial-mesenchymal transmission (EMT; P value<0.005).
As there was a significant enrichment in the immune/ inflammatory pathway through GSEA analysis, we conducted immune infiltration analysis. ESTIMATE algorithm displayed significant differences in the immune score (P=0.02) and ESTIMATE score (P=0.027) between high and low autophagy risk groups in the TCGA CRC cohort (Figure 7A). Infiltration of plasma cells and Tregs were enriched significantly in the high autophagy risk group compared with the low autophagy risk group in the GSE39582 and TCGA cohort (Figure 7C). By contrast, M1 macrophage infiltration turned out to be significantly lower in the high autophagy risk group (Figure 7C).
3.4 Adjuvant chemotherapy effects on different autophagy risk groups
Survival analysis conducted for different autophagy risk groups with and without adjuvant chemotherapy respectively, showed that for early stage CRC patients without adjuvant chemotherapy in the high autophagy risk group displayed poorer DFS in the GSE39582 cohort (Figure 8A, 5[2.38–10.5]; P=2.42×10-6). However, DFS showed no significant difference between high and low autophagy risk groups for early CRC patients with adjuvant chemotherapy (Figure 8B, P=4.46×10-1). Similar results were observed in the TCGA cohort (Figure 8C, 2.05[0.9–4.67]; P=8.22×10-2 and Figure 8D, P=5.84×10-1).