2.1 Differentially expressed IRGs

All original data were obtained from the TCGA database. We examined the differentially expressed genes between PRAD and corresponding normal tissues to screen out IRGs involved in oncogenesis. The abnormally expressed genes were identified using the software edgeR (version 4.0.0, John Chambers and colleagues, Auckland, New Zealand) (https://www.r-project.org/), as described previously (adjusted p value < 0.05 and |log2 (fold change)| > 1). The intersection between differentially expressed genes and the IRGs was termed differentially expressed IRGs.

2.2 Pathway and functional enrichment analyses

We carried out Pathway and functional enrichment analyses to investigate the biological significance of these differentially expressed IRGs. Kyoto Encyclopedia of Genes and Genomes (KEGG) [6] and Gene Ontology (GO) [7] enrichment analyses were carried out. Terms and pathways were regarded as significantly enriched objects if their false discovery rate was less than 0.05. The package ggplot2 was used to carry out the visualization. Pathway illustrations were generated using the “pathview” R package.

2.3 Constructing an IRGs-related immune signature for PRAD

We identified immune genes associated with prognosis. To verify their robustness and predictive ability, we developed a prognostic risk model. We first screened for immune-related genes with prognostic significance by univariate Cox proportional hazards regression analysis. Genes with a p-value of less than 0.05 were examined by the glmnet software using the minimum absolute contraction and selection operator (Lasso) to prevent over-fitting [8]. After screening by the Lasso model, a multivariate Cox proportional hazards model was employed to develop the immune-related risk model as follows: risk score = Gene A level × coefficient a + gene B level × coefficient B + gene C level × coefficient C + … + gene N level × coefficient N. The prognosis for PRAD patients is represented by the risk score in the model, where the lower the risk score, the better the prognosis. The patients were classified into high-risk and low-risk groups using the median risk score as the cut-off value. With log-rank p-value < 0.05 considered statistically significant (for packet survival and investigator), predictive ability was calculated using Kaplan-Meier survival curves. Using a time-dependent receiver operating characteristic curve (ROC), the prediction ability of this immune signature was evaluated [2, 9, 10].

2.4 Calculating the ratio of tumor-infiltrating immune cells

By using a deconvolution technique, CIBERSORT could determine the ratios of infiltrated immune cells from tissue transcriptional profiles. We computed the ratios of 22 types of tumor-infiltrating immune cells using the TCGA-PRAD transcriptional profiles and the CIBERSORT R script [11, 12].

2.5 Statistical analysis

Chi-square and Student’s t-tests were conducted to examine differences among variables. The effects of multiple clinic-pathological characteristics along with the immune signature on patients’ survival were evaluated using univariate and multivariate cox regression analyses [13]. Using R software programming language (version 0.0.0). Statistical analysis was carried out. The package “pheatmap” was used to create the heat maps. Statistical significance was defined as p-value < 0.05.

2.6 Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

The specific experimental steps are as described in the previous paper [14]. In brief, total RNA was extracted from normal tissue and tumor tissue from PRAD patients using TRIzol (15596026, Thermo Fisher Scientific, Waltham, MA, USA). Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was conducted on the obtained RNA from each sample (2 μg) with FastStart Universal SYBR® Green Master (4913914001, Roche, Basel, Switzerland., USA) on a Q5 PCR System (Roche, USA). The primers were showed in Table 1.

3.1 Identification of prostate cancer prognosis-related differentially expressed genes

We first classified the prostate cancer samples in the TCGA database as “Alive” and “Dead” according to patients’ overall survival profiles. After adjusting the cut-off value adjusted p value < 0.05 and |log2 (fold change)| > 1). Then, 1309 Differentially expressed genes were screened. At the same time, 2489 immune-related genes with an “Immune” correlation score greater than 2 were screened from the GeneCards database. After crossing, 18 immunologically related prognostic differentially expressed genes were obtained.

3.2 The prognosis of PRAD patients can be predicted by the five-gene immune signature

To evaluate the predictive significance of the chosen 18 differentially expressed IRGs and to avoid over-fitting, we further performed a LASSO analysis in which five of the 18 genes were predictors of patient outcome (Fig. 1A,B). We ranked patients’ risk scores and outlined their distribution (Fig. 1C Top). A dot plot (Fig. 1C Middle) was used to display the risk scores, survival status and overall survival of the patients. The heat map displays the five IRGs’ expression patterns in patients with varying risk scores. In the low-risk group, Insulin Like Growth Factor 1 (IGF1) and Opioid Receptor Kappa 1 (OPRK1) were expressed at high levels, whereas in the high-risk group, Anti-Mullerian Hormone (AMH), Defensin Alpha 5 (DEFA5) and Serum Amyloid A1 (SAA1) were expressed at high levels. The predictive ability of this immune signature was further verified. Prostate cancer patients were classified into high-risk and low-risk groups (categories). Those with low-risk scores had longer Overall survival (OS) duration, according to the survival analysis (p = 0.004) (Fig. 1D). Moreover, this immune signature’s area under the curve (AUC) value was 0.875 (Confidence interval (CI): 0.796–0.953) (Fig. 1E). Time-dependent ROC curves at 1-year, 3-year and 5-year durations (AUC = 1.000, 0.891 and 0.912, respectively) are also displayed (Fig. 1F).

Fig. 1.Lasso analysis and forest plot show multivariate Cox model results for 5 immune-related genes. (A) Cross validation of optimal parameter selection in LASSO model. (B) LASSO coefficient profiles of 18 prognostic immune-related genes; (C) The risk factor graph shows the distribution of risk scores (Top), the relationship between risk scores and survival time (middle), and the expression patterns of five immune-related genes in the high-and low-risk groups (bottom); (D) Kaplan-Meier curves of overall survival of PRAD patients in high- and low-risk score groups; (E) Receiver performance curve based on risk score; (F) Time-dependent receiver operating characteristic curves. AUC: area under the curve; TPR: True positive rate; FPR: False positive rate; OS: Overall survival; AMH: Anti-Mullerian Hormone; DEFA: Defensin Alpha; IGF: Insulin Like Growth Factor; OPRK: Opioid Receptor Kappa; SAA: Serum Amyloid A.

3.3 The tumor immune microenvironment and the 5-IRG signature

To generate histograms displaying the proportion of 22 immune cells infiltrating in tissue samples from the PRAD patients, we computed the proportion of 22 immune cell types in each PRAD sample using the CIBERSORT method (Fig. 2A,B). The proportion of immune cells that differed between the groups with high and low-risk scores was compared. The percentage of B cells, Eosinophils, Macrophages, mast cells, Neutrophils, NK (Natural killing cell) CD56 (Neural Cell Adhesion Molecule 1) bright cells, NK cells Tgd, PDC (plasmacytoid dendritic cell), Th1 cells and T helper cells differed remarkably (Fig. 2C) [15].

Fig. 2.Evaluation of immune infiltration. (A,B) Stacked bar graphs showing the proportion of infiltrating immune cells between the high-and low-risk groups. (C) Boxplot of the difference in immune cell infiltration between the high-and low-risk groups (Wilcoxon’s test; *: p < 0.05; **: p < 0.01; ***: p < 0.001). Horizontal coordinates indicate the type of infiltrating immune cells, and the vertical axis indicates the proportion of infiltrating immune cells. Immune cell abbreviation: Th: Helper T lymphocytes; TR: regulatory T cell; NK: Natural killing cell; CD56: Neural Cell Adhesion Molecule 1; DC: Dendritic cell; PDC: plasmacytoid dendritic cell; TFH: Follicular helper T cells; Tgd: gamma delta T.

We compared the expression differences of prognostic DEirgs in TCGA-PRAC, with only AMH being significantly overexpressed in the tumor (Fig. 3A) while having the same expression differences in paired samples of tumor tissue (Fig. 3B). The expression of the five genes was significantly different and were demonstrated using a between-group differential gene volcano map (Fig. 3C) and a sequence map (Fig. 3D) based on patient survival status. Subsequently, both univariate and multivariate Cox regression analyses were carried out to investigate the prognostic value of model DEIRGs (Fig. 3E,F & Table 2).

Fig. 3.Differential expression of prognostic DEIRGs in TCGA-PRAC. (A,B) The same difference was found in tumor tissues and normal samples. (C,D) Volcano plot and sequence plot showed the difference and significance of the five genes in the prognostic model based on the survival status of the patients. (E,F) Forest plot of 5 modifiers with prognostic value in univariate and multivariate Cox regression models. DEIRGs: Different expressed immune-related genes; TCGA: The Cancer Genome Atlas; PRAC: Prostate adenocarcinoma; AMH: Anti-Mullerian Hormone; DEFA5: Defensin Alpha 5; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1; GAPDH: Glyceraldehyde-3-Phosphate Dehydrogenase.

3.4 Relationship between survival prognosis and 5-gene expression difference in prostate cancer patients

We explored the association of differences in 5-gene expression and disease-specific survival, overall survival and progression-free interval in prostate cancer patients. Although high AMH expression was not significantly linked to low OS (overall survival) rate. Prostate cancer patients with low expression of SAA1 had poor Overall Survival (p = 0.041) (Fig. 4A), patients with elevated AMH expression levels had worse overall survival (p = 0.127) (Fig. 4B). And patients with OPRK1-low expression had poor disease-specific survival and overall survival (p = 0.033 and 0.022, respectively) (Fig. 4C,D); Progression-free interval and disease-specific survival (p = 0.004, 0.029, respectively) were obtained in patients with prostate cancer with low expression of IGF1 (Fig. 4E,F) [16].

Fig. 4.The relationship between overall survival rate and gene expression of SAA1 (A), AMH (B), OPRK1 (D). The relationship between disease-specific survival rate and gene expression of OPRK1 (C), IGF1 (F). The relationship between progression-free interval and gene expression of IGF1 (E). AMH: Anti-Mullerian Hormone; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1.

3.5 The efficacy of prognostic identification was verified

To better understand how a model gene’s expression affects clinical prognosis in prostate cancer, ROC curve analysis was used to verify its efficacy in differentiating different clinical outcomes. We first validated the diagnostic efficacy of a five-gene prognostic model for prostate cancer disease status as DEFA5 (AUC = 0.464), IGF1 (AUC = 0.636), SAA1 (AUC = 0.553), OPRK1 (AUC = 0.493), and AMH (AUC = 0.809), respectively (Fig. 5A). And then, we validated the efficacy of 5-gene prognostic models for the identification of Clinical variable outcomes of Alcohol history, Smoker, Lymphnode neck dissection, Radiation therapy, Age, Clinical stage, Histologic grade, Clinical T stage, pathologic stage, pathologic N stage, pathologic T stage, respectively (Fig. 5B–L). These results suggest that, for the assessment and anticipation of clinicopathological factors in prostate cancer, the 5-gene prognostic model may be a potential biomarker.

Fig. 5.ROC curves validate diagnostic efficacy. (A) The diagnostic efficacy of the 5-gene prognostic model for prostate cancer; (B–L) To verify the efficacy of the 5-gene prognostic model in the identification of Clinical variable outcomes of Alcohol history, Smoker, Lymphnode neck dissection, Radiation therapy, Age, Clinical stage, Histologic grade, Clinical t stage, Pathologic stage, Pathologic n stage, Pathologic T stage, respectively. AUC: area under the curve; AMH: Anti-Mullerian Hormone; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1. DEFA5: Defensin Alpha 5; TPR: True positive rate; FPR: False positive rate.

3.6 The prognostic model is associated with prostate cancer

To further explore the association of 5-gene prognostic models with prostate cancer, we explored their association with the expression of well-known molecular markers known to be associated with prostate cancer. We selected some reliable prostate cancer biomarkers for further analysis to explore whether the 5-gene prognostic model is associated with the molecular mechanisms of prostate cancer initiation and progression. It has been reported that transmembrane serine endopeptidase 2 (TMPRSS2) and ETS-related genes (ERG) frequently fuse during Principal component analysis (PCA) [17]. A mitochondrial and peroxisome enzyme, Alpha-methylacyl-coa racemase (AMACR), participates in fatty acid oxidation and intermediate products of bile acids [18]. The over-expression of AMACR on prostatic epithelial cells occurs in about 80% of prostate cancer cases, and the detection of AMACR by immunohistochemistry (IHC) helps to distinguish benign prostate cancer from prostate cancer [19]. At the same time, studies have reported that AMACR overexpression in prostate cells predicts poor prognosis (associated with a high Gleason score and high initial PSA levels) and an increased risk of bone metastasis. In addition, there have been reports linking several Solute Carrier Family 45 Member (SLC45A) family genes to prostate cancer [20, 21, 22]. As SLC45A3 (prostate cancer-associated protein 6/P501S/protein) is mainly expressed in the Gleason apparatus of prostate epithelial cells and is highly specific for prostate gland cells, it has been used to distinguish between other tumor types and metastatic prostate cancer [22, 23]. SLC45A3 may be expressed in PSA-negative prostate tumors, so the combined use of these markers may increase the sensitivity to identify prostate cancer metastases. SLC45A3 is weakly expressed in some aggressive tumors and is associated with increased Gleason score and risk of tumor recurrence. Overall findings revealed that in the 5-gene prognostic model, OPRK1 expression was positively linked to BRCA1 DNA Repair Associated (BRCA1), BRCA2 DNA Repair Associated (BRCA2), and AMACR (Fig. 6A,B,E). The positive correlation between OPRK1 expression and TMPRSS2 was significant, while the negative correlation between SAA1, AMH and TMPRSS2 was significant (Fig. 6C). SAA1 was positively correlated with Hexokinase 3 (HK3), SLC45A1, SLC45A3, SLC45A4 (Fig. 6D,F–H).

Fig. 6.Prognostic model is associated with prostate cancer. The overall results showed that the expression of OPRK1 in five prognostic model molecules was correlated with the expression of BRCA1, BRCA2 (A,B), TMPRSS2 (C), HK3 (D), AMACR (E), SLC45A1, SLC45A3, SLC45A4 (F–H) in lollipop graph, the size and color of the circle represent significance. BRCA1: BRCA1 DNA Repair Associated; AMH: Anti-Mullerian Hormone; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1; DEFA5: Defensin Alpha 5; HK3: Hexokinase 3; TMPRSS2: transmembrane serine endopeptidase 2; AMACR: Alpha-methylacyl-coa racemase.

3.7 Functional enrichment analysis

In addition, we summarised the differential genes between the low-risk and high-risk categories, followed by GO and KEGG enrichment analyses. Based on the results of the GO enrichment analysis (Fig. 7A–D & Table 3), the BP with significance were negative regulation of the reproductive process, ovulation cycle, oogenesis, interleukin-1 production and regulation of Interleukin-1 production, respectively. Cell component (CC) was enriched in transport vesicles, exocytic vesicles, cytoplasmic vesicle lumen, vesicle lumen and secretory granule lumen. Hormone activity, growth factor activity, insulin receptor binding, insulin-like growth factor receptor binding, and transforming growth factor beta receptor binding were all enriched in MF. Furthermore, KEGG analysis revealed some immune-related pathways.

Subsequently, the 5 genes of the prognostic model were mapped in red letters on the KEGG pathway map. Pathway enrichment analysis suggested that the 5-gene prognostic model was enriched in aldosterone regulation, steroid production, long-term depression, longevity regulatory pathways, and Staphylococcus aureus infection, and the well-known Epidermal Growth Factor Receptor (EGFR), Hypoxia Inducible Factor 1 (HIF-1), The AMP-activated protein kinase (AMPK), and other signaling pathways, suggesting that a variety of cancer-related (Fig. 7E,F & Table 4). Fig. 8 shows the role of 5 genes involved in Prostate cancer, Long-term depression, Endocrine resistance, Transforming Growth Factor Beta (TGF-β) signaling pathway, and EGFR tyrosine kinase inhibitor resistance.

Fig. 7.GO and KEGG enrichment analysis of 5-gene prognostic model molecules. (A–D) GO enrichment analysis bar graph, bubble graph, chord graph and cycle graph; (E,F) Kobas database was used to analyze the KEGG pathway enrichment results of the 5-gene prognostic model. GO: Gene Ontology; BP: Biological process; CC: Cell component; MF: Molecular function; AMH: Anti-Mullerian Hormone; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1; DEFA5: Defensin Alpha 5; TGF-β: Transforming Growth Factor Beta; cAMP: Cyclic adenosine monophosphate.

Fig. 8.KOBAS database was used to analyze the KEGG pathway of 5-gene prognostic model. (A) Prostate cancer, (B) Long-term depression, (C) Endocrine resistance, (D) Epidermal Growth Factor Receptor (EGFR) tyrosine kinase inhibitor resistance, (E) Transforming Growth Factor Beta (TGF-β) signaling pathway, (F) p53 signaling pathway.

3.8 Gene set enrichment analysis based on a 5-gene prognostic model to predict differences between low and high-risk groups for prostate cancer

The potential mechanisms by which prognostic models regulate OS timing in PRAD patients need to be further explored. To find the possible signaling pathways involved in the prognostic model in the two risk groups, we performed GSEA (Fig. 9, Supplementary Table 1). The KEGG high-risk group was primarily enriched in Ribosome, Cell Cycle, Spliceosome, Ascorbate and Aldarate Metabolism, and Base Excision Repair, whereas the low-risk group was primarily enriched in Cardiac Muscle Contraction, Hypertrophic Cardiomyopathy HCM, Dilated Cardiomyopathy, Proximal Tubule Bicarbonate Reclamation, and Primary Immunodeficiency (Fig. 9A,B). In WikiPathways, high-risk groups were enriched in Striated Muscle Contraction Pathway, and Extrafollicular B Cell Activation by Sarscov2, Folate Metabolism, Selenium Micronutrient Network, and Complement System (Fig. 9C,D). The high-risk groups in the Reactome database were mostly enriched for striated muscle contraction, CD22-mediated B-cell receptor (BCR) regulation, clearance of heme in plasma, and creation of initial triggers for the complement system, such as C4 and C2 activators; The low-risk group was enriched in Striated Muscle Contraction, CD22 Mediated BCR Regulation, Scavenging of Heme from Plasma, spindle checkpoint and sister chromatid, of C 4 and C 2 creation Activators, Initial Triggering of Complement (Fig. 9E,F). In Persistent Identifier (PID), the high-risk group was enriched in Polo Like Kinase 1 (PLK1), ATR Serine/Threonine Kinase (ATR), Aurora B, Fanconi and Aurora a Pathway, while the low-risk group was enriched in NCADHERIN, INTEGRIN3, INTEGRIN2, Interleukin (IL42) pathways (Fig. 9G,H). Finally, the expression of hub gene was confirmed by PCR (Fig. 10).

Fig. 9.Gene set enrichment analysis based on a 5-gene prognostic model to predict differences between high-and low-risk groups for prostate cancer. Significant enrichment entries were included in the REACTOME (A,B), WP (C), PID (D), Kegg (E) and BIOCARTA (F) databases, respectively. KEGG: Kyoto Encyclopedia of Genes and Genomes; PID: Persistent Identifier. ATR: Serine/Threonine Kinase; CD56: Neural Cell Adhesion Molecule 1; PID: Persistent Identifier; PLK1: Polo Like Kinase 1; IL42: Interleukin; BCR: B-cell receptor.

Fig. 10.The gene expression of hub genes by PCR. AMH: Anti-Mullerian Hormone; DEFA5: Defensin Alpha 5; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1.