3.1 Data processing and identification of DEARGs
After the obtained data were subjected to background correction, data normalization, and removal of the difference between batches, a total of 222 autophagy-related genes and 120 DEARGs were obtained.
3.2 DEARGs function and pathway enrichment analysis
The results of GO analysis demonstrated that the biological process (BP) changes of the DEARGs were significantly enriched in ‘autophagy’, ‘process utilizing autophagic mechanism’, ‘macroautophagy’, ‘autophagosome assembly’, ‘autophagosome organization’, ‘regulation of autophagy’, ‘vacuole organization’, and ‘autophagy of mitochondrion’; Cell components (CC) were mainly enriched in ‘autophagosome’, ‘phagophore assembly site’, ‘autophagosome membrane’, ‘vacuolar membrane’, ‘melanosome’, ‘pigment granule’, ‘phagophore assembly site membrane’, and ‘late endosome’; Molecular function (MF) mainly involved ‘ubiquitin protein ligase binding’, ‘ubiquitin−like protein ligase binding’, ‘chaperone binding’, ‘protein serine/threonine kinase activity’, and ‘heat shock protein binding’. KEGG pathway enrichment was mainly enriched in ‘autophagy’, ‘shigellosis’, ‘kaposi sarcoma-associated herpesvirus infection’, ‘human cytomegalovirus infection’, ‘mitophagy’, ‘protein processing in endoplasmic reticulum’, and ‘apoptosis’ (Figure 1A).
3.3 PPI network construction, hub autophagy-related gene selection, and functional similarity analysis of DEARGs
The results of the PPI network showed that three genes, including TP53, ATG16, and MTOR, had larger weights and stronger correlations (Figure 2A). The interaction maps of the 10 hub autophagy-related genes (BECN1, ATG7, PIK3C3, ATG12, GABARAPL1, GABARAPL2, SQSTM1, MAP1LC3A, MAP1LC3B, and GABARAP) were screened (Figure 2B and Table 1).
To further identify the strength of the interactions between the hub autophagy-related genes, we sorted them based on their average functional similarity, and the results showed that MAP1LC3B, GABARAPL2, GABARAPL1, and GABARAP were the genes with the closest interactions. Among them, MAP1LC3B is the only gene with a cut-off value > 0.7, with the highest average functional similarity (Figure 2C).
3.4 Construction of prognostic autophagy-related genes signature
To study the prognostic value of DEARGs, the univariate Cox analysis revealed that 25 genes including BNIP3, MYC, VEGFA, ATG9B, NKX2.3, ATG4A, BAG1, CALCOCO2, ATF4, ATG7, ATG4D, AMBRA1, BAK1, APOL1, EGFR, CCL2, MAPK1, CASP1, MAP2K7, DLC1, DNAJB9, VAMP3, PEX3, CTSD, and BCL2 had a prognostic value for OS (Figure 3 and Supplementary material 1). To further improve the accuracy of the results, we performed Lasso Cox regression analysis on the 25 DEARGs. The results showed that 9 genes (BNIP3, MYC, BAG1, CALCOCO2, ATF4, AMBRA1, EGFR, MAPK1, and PEX3) were correlated to prognosis (Figure 4 and Table 2). These 9 genes were subsequently used to construct a gene signature. We calculated the risk score of each sample, and divided the samples into high-risk and low-risk groups using the median risk as a cut-off value, and obtained a risk factor correlation diagram of the prognostic autophagy-related gene risk score. Compared to the low-risk group, the number of deaths in the high-risk group was significantly greater, and the levels of ATF4, BAG1, MYC, and BNIP3 were higher in the high-risk group, while the expression levels of MAPK1, EGFR, CALCOCO2, and AMBRA1 were lower (Figure 5A). Kaplan-Meier survival analysis showed that the overall survival rate of patients in the high-risk group was significantly lower than that in the low-risk group (P <0.0001) (Figure 5B).
3.5 Relationship between prognostic autophagy-related gene signature and clinical characteristics
As shown in Figure 6 and Supplementary material 2, the risk score was significantly associated with the distant metastasis of OS (P = 0.05), but not with age (P = 0.710) and gender (P = 0.887).