Differentially expressed immune-related genes
The research process was showed in Figure 1. Altogether RNA-seq and clinical data of 499 PCa tissue samples and 52 non-tumor samples were downloaded from TCGA. Among these patients, a total of 499 primary PCa patients with gene expression data and clinical follow-up information was involved in the current study. According to the result, we initially screened 193 differentially expressed genes (Figure 2A-2B).
Network analysis of TF-IRG interaction
To explore the interaction between transcription factors (TF) and IRGs, we extracted 22 differentially expressed TFs by analyzing the RNA-seq data described above. The heatmap of TFs with markedly different expression between the PCa tissues and normal prostate tissues were shown in Figure 3A, and we obtained 11 up-expressed and 11 down-expressed TFs (Figure 3B-C). The TF-IRGs interaction network was detailed in Figure 3D.
Functional annotation of the differentiallyexpressed immune-related genes
Functional enrichment analysis of the 193 differentially expressed IRGs offered that the biological understanding of these genes. The GO terms function and KEGG pathway enrichment of these genes were summarized in Table 1. According to the results of DAVID, we found that the top enriched GO terms for biological processes were: positive regulation of response to external stimulus, positive regulation of cell migration and leukocyte migration, and for cellular components were:immunoglobulin complex, circulating, external side of plasma membrane, and immunoglobulin complex. On the basis of molecular function, genes were mostly enriched in terms of receptor ligand activity, growth factor activity, and cytokine activity. The overview schematic of the analysis results is displayed in Figure 4. Besides, in the KEGG pathway enrichment analysis for the differentially expressed immune-related genes, these genes were shown to be notably associated with Pathways in Cytokine-cytokine receptor interaction, Ras signaling pathway, Neuroactive ligand-receptor interaction, MAPK signaling pathway, EGFR tyrosine kinase inhibitor resistance, and so on. As shown in Figure 4A, the Z-score of enriched pathways less than zero indicated that most of the cancer pathways were more likely to be decreased (Figure 4B).
Construction and validation of the immune-related risk signature
These survival-related genes were subjected to a univariate Cox regression analysis and multivariate Cox regression analysis to remove the genes that might not be an independent indicator in prognosis monitoring. Finally, several prognostic IRGs were obtained. The relationships between the expression profiles of 193 differentially expressed IRGs were obtained from TCGA, resulting in six prognosis-related IRGs(Figure 5A). In order to improve the robustness, six prognosis-related immune-related (S100A2, NOX1, IGHV7-81, AMH and AGTR1, BIRC5) were selected for further multivariate Cox regression model by SPSS 24.0 (Figure 5B). However, the gene of AGTR1 showed no significant prognostic value with P>0.05. Finally, five genes including S100A2, NOX1, IGHV7-81, AGTR1 and AMH were identified to develop the model (Table 2).
Evaluation of the immune-related prognostic signature for prostate cancer
After constructing our prognostic model, we validated and evaluated this model. The results from K-M analysis indicated that the high-risk group had shorter survival time than the low-risk group, and p value of the log-rank test was 1.163e−03, and the heatmap of prognosis-related IRGs were also displayed(Figure 6A-6B). To evaluate how well the immune-related prognostic signature predicts the prognoses of PCa, the time-dependent ROC curve analysis was carried out. The receiver operating characteristic curves (ROC) were created by the“survivalROC” package, and the area under the curve (AUC) values was calculated to evaluate the specificity and sensitivity of the model. The AUC for the prognostic signature was 0.985, demonstrating the competitive performance of the prognostic signature for survival prediction in the TCGA dataset(Figure 6C). The riskscore distribution of patients in different risk groups showed that the high-risk group has high risk score, and survivalstate program showed that patients with high riskscore have high dead rates(Figure 6E-6F). A prognostic nomogram was also performed to visualize the relationship between nomogram and survival rates in patients with immune-related genes based on the Cox proportional hazard regression model by means of “rms” package of R software(Figure 6D).
Integrated prognostic signature by combining with clinical parameters
Results showed that high S100A2 gene expression has a lower biochemical recurrence, high AMH gene expression has a higher gleason score, lymph node metastasis rate and tumor grade, a lower lymphnodes positive, high ATGR1 gene expression has a lower psa value(Figure 7A-7F).