Identification of an autophagy-related risk signature for the prognosis of prostate cancer
The prognostic value of ARGs was performed by univariate COX regression in 499 prostate tumor samples in the TCGA database. 17 genes (ATG16L1, FADD, GABARAPL2, NKX2-3, MYC, MAPK8IP1, WDR45B, MTMR14, HGS, USP10, NPC1, BIRC5, BNIP3, ATG3, RAB24, ULK3 and RUBCN) were identified to be significantly correlated with OS in prostate tumor samples. Among the 17 genes, all the genes were identified as risk factors (p < 0.05, HR > 1, Figure 1A). And then a multivariate Cox regression was performed to develop the following autophagy-related risk signature related to the survival of PCa patients. The calculation formula of risk score is as follows [32]: Risk score = (FADD expression×2.341572) + (GABARAPL2 expression×4.34965) + (MYC expression×2.430837) + (RAB24 expression×3.570068) + (RUBCN expression×4.857469) + (NPC1 expression×1.739033). Then patients were divided into low-risk group (n=177) and high-risk group (n=176) using the median risk score as a cutoff value. Our data showed that the survival time of high risk group is significantly shorter than the low risk group (Figure 1B). Risk genes showed significant expression patterns according to the risk value (Figure 1C) and all of the deaths occurred in the high-risk group (Figure 1D). The heatmap of the 6 ARGs expression levels in the TCGA dataset, high-risk group expressed higher levels of risky genes (Figure E). ROC curves for 3-year overall survival were performed to evaluate the predictive power of the six-gene risk signature (Figure 1F). The 3-year AUC of our signature was 0.928, which was obviously higher than that of age (AUC = 0.580), tumor T stage (AUC = 0.641), and tumor N stage (AUC = 0.694).
The correlation between the autophagy-related risk signature and clinical factors in prostate cancer
An analysis was applied to compare the correlation between predicted 6-gene signature and conventional clinical factors of prostate cancer, Results showed that MYC expression was correlated with patients’ age (P = 0.009). FADD expression was correlated with both T classification (P = 0.009) and lymphatic invasion (P = 0.012). NPC1 expression was associated with lymphatic invasion (P = 0.017). RUBCN expression was correlated with lymphatic invasion (P = 2.098e-04). We further analyzed the correlation between the risk score and these clinical factors, we found that T classification was correlated with the risk score (P = 0.024) (Figure 2).
Identification of DEMs and DEGs in prostate cancer
There are 21 pairs of PCa samples and matched adjacent non-tumor prostate samples were collected and processed for microRNA detection, and mRNA expression analysis was performed by using 30 matched malignant and normal prostate tissue samples from 15 prostate cancer patients. All the samples were obtained from GEO. Results showed that a total of 12 DEMs and 1073 DEGs were detected. Considered as the criteria of |log FC| >1 and adjust P value < 0.05, we finally identified 2 up-regulated and 10 down-regulated DEMs (Figure 3A and C). Meanwhile, 413 up-regulated and 660 down-regulated DEGs were extracted (Figure 3B and D).
Enrichment analysis of the DEMs and DEGs
To improve our understanding of the biological information of the 12 DEMs in prostate cancer, we performed GO annotation and biological pathway analyses by using the software of FunRich. Regarding BP, the DEmiRNAs were significantly enriched in regulation of nucleic acid metabolism, , nucleotide, nucleoside, and nucleobase (Figure 4A). Regarding CC, the DEmiRNAs were significantly enriched in nucleus, cytoplasm, lysosome, golgi apparatus, and endocytic vesicle membrane (Figure 4B). In addition, the significantly enriched GO terms in MF was transcription factor activity (Figure 4C). As shown in Figure 4D, the pathways of biological processes were TRAIL signaling pathway, Class I PI3K signaling events mediated by Akt, PDGFR-beta pathway, mTOR signaling pathway, EGF receptor pathway, VEGF and VEGFR network, IFN-gamma pathway, ErbB receptor signaling network, Glypican pathway, and PDGF receptor signaling network. Alike, all 1073 DEGs were also uploaded to the FunRich, the results of GO analysis indicated that 1) for BP, DEGs were significantly enriched in cell communication, signal transduction, and cell growth and/or maintenance; 2) for CC, DEGs were particularly enriched in the basement membrane, extracellular region, proteinaceous extracellular matrix, plasma membrane, extracellular space, extracellular matrix, and extracellular. 3) for MF, DEGs was only enriched in extracellular matrix structural constituent (Figure 4E-4G). Additionally, biological pathway analysis showed these DEGs were mostly enriched in mesenchymal-to-epithelial transition and epithelial-to-mesenchymal transition (Figure 4H).
miRNA–mRNA network
FunRich software was utilized to predict potential target genes from DEMs. All the 12 DEMs were inducted into the FunRich software. There were 1980 target genes found. Then, we further assessed the intersection of 1980 target genes and 1073 DEGs, and obtained 104 overlapping genes for subsequent analysis (Figure 5A and B). The network results of DEMs and overlapping genes were calculated using FunRich software, and the results were visualized in Cytoscape software (Figure 5C). Notably, hsa-miR-148a targeted 22 genes, including ADAMTS18, ADAMTS5, CAV2, CCDC85A, COL6A3, DNAJB4, EMX2, FBN1, FOXF1, GPM6A, HLF, MYBL1, NDP, PRICKLE2, S1PR1, SULF1, TSPAN18, ZNF804A, B4GALT6, COL4A1, LAMA4, TGFB2; hsa-miR-133b targeted 3 genes, including SH3GL2, SFXN2, CDCA8; hsa-miR-204 targeted 5 genes, including SLC43A1,SGIP1, PRR15L, SFXN2, EPHA5; hsa-miR-222 targeted 2 genes, including STMN1 and SBK1; hsa-miR-221 targeted 2 genes, including STMN1 and SBK1; hsa-miR-31 targeted 4 genes, including MBOAT2, PPP1R9A, CTNND2, PRSS8; hsa-miR-205 targeted 5 genes, including DSC2, NKX2-3, HS3ST1, ACSL1, EPB41L4B; hsa-miR-455 targeted 5 genes, including COL2A1, HOXC4, COLEC12, KLK12, STEAP2; hsa-miR-145 targeted 5 genes, including HOMER2, IGSF5, LDLRAD3, PGM3, TMEM178A; hsa-miR-375 targeted 2 genes, including ISL2 and KCNE3; hsa-miR-376c targeted 2 genes, including ALCAM and NKX3-1; These results indicated that in all the relationships between DEMs and DEGs, miR-205 can specifically regulate the expression of autophagy-related gene NKX2-3.
Construction of a interaction network between NKX2-3 and autophagy-related risk signature
STRING database was performed to diagram a network of interacting relationships among the 17 ARGs obtained in Figuge 1A. As shown in Figure 6A, ATG16L1 can act as a bridge node between NKX2-3 and three of the six genes risk signature (GABARAPL2, RUBCN, RAB24). In other words, we speculated whether NKX2-3 can affect above risk signature through regulating GABARAPL2, RUBCN, and RAB24 expression. To validate this network further, Pearson correlation analysis was used to determine the correlation between NKX2-3 expression and the three genes associated with the risk signature. The results indicated that significant positive correlations were observed between NKX2-3 and ATG16L1, GABARAPL2, KIAA0226 (RUBCN) respectively (Figure 6B-E).
Identification of differentially expressed ARGs in prostate cancer
Altogether clinical data and transcriptome expression profiles were downloaded from TCGA. Expression values of 232 ARGs were extracted, we finally identified 5 up-regulated and 8 down-regulated ARGs under cut-off criteria of FDR<0.05 and |logFC|>1 (Figure 7A and B). Furthermore, as presented in Figure 7C, box plots displayed expression patterns of 8 down-regulated genes (NRG2, BCL2, NRG1, HSPB8, FAM215A, TMEM74, TP63, and ITPR1) and 5 up-regulated genes (NKX2-3, CDKN2A, BIRC5, CAMKK2, and ATG9B).
Survival analysis of NKX2-3 and clinical parameters in patients with prostate cancer
Based on the results described above, we found that the abnormally high expression of autophagy-related gene NKX2-3 in prostate cancer not only serves as a prognostic risk factor for prostate cancer patients, but also has a specific regulatory relationship with miR-205. Therefore, we further analyzed the clinical value of NKX2-3. The cBioPortal is an open-access resource for interactive exploration of multidimensional cancer genomics data sets and provides comprehensive analyses [33, 34]. The NKX2-3 overall survival analysis using cBioPortal datasets termed prostate cancer (TCGA, Firehose Legacy). Our data showed that the altered group had a shorter overall survival time than the unaltered group (Figure 8A). Then, we further analyzed the relationship between NKX2-3 and some common clinical parameters. The results showed that the lymph nodes examined number and distant metastasis rate were significantly increased in altered group, compared with the corresponding unaltered group (Figure 8B and C).
GSEA by NKX2-3 stratification in prostate cancer
To explore the mechanism of NKX2-3, patients were separated into high/low expression groups based on the median expression of NKX2-3 and then were subjected to GSEA. The GSEA analysis revealed that the high-expression group of NKX2-3 has strikingly upregulated genes enriched in Myc Targets V1, Unfolded Protein Response and Myc Targets V2 (Figure 9).
NKX2-3 may serve as a potential predictor for the efficacy of anti-PD-1 therapy in prostate cancer
To further test the correlation of NKX2-3 with immunotherapy response of prostate cancer, we further analyzed the relationship between NKX2-3 and PD-1 (programmed cell death 1, PDCD1). The relevance between the expression of NKX2-3 and PD-1 was analyzed by TIMER, which is a comprehensive resource for systematical analysis of immune infiltrates across diverse cancer types [35, 36]. Our data showed that as the expression of NKX2-3 increased, the expression of PD-1 decreased, thus NKX2-3 was negatively correlated with PD-1 in prostate cancer (Figure 10A). Furthermore, as shown in Figure 10B, a positive correlation was detected between NKX2-3 expression and the TMB (P<0.05; (tumor mutation burden) of prostate cancer (P = 1.3e-06).