The transcriptome profiling and clinical data of PCa were obtained from TCGA (https://portal.gdc.cancer.gov) with TCGAbiolinks (HTSeq-Counts). Two other datasets (GSE21034 and GSE70770) were accessed from GEO (https://www.ncbi.nlm.nih.gov/geo/). The TCGA-PRAD dataset was selected by the following steps: (1) The follow-up was more than 20 days. (2) Samples without complete BCR follow-up clinical information were removed. We obtained 411 PCa patients with complete BCR follow-up information and 406 patients with complete clinical information ( Table S1 ). The EdgeR package was used to analyze differentially expressed genes (DEGs) between tumor and paracancerous tissue on the R 4.1.3 platform, and FDR < 0.05 and |log2-fold change| ≥ 1 were considered to be statistically significant. ERS-associated genes were collected from GeneCards (https://www.genecards.org/), and the correlation score ≥ 7 was selected, as Huang et al. reported (12).

Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analysis (GSEA) were used to perform enrichment analysis of the DEGs with the “ClusterProfiler” package (13). Gene Set Variation Analysis was carried out with “GSVA” packages.

The ConsensusClusterPlus package (14) was used to perform unsupervised hierarchical clustering to identify differentially expressed ERS-related clusters via pam algorithms. A total of 1,000 iterations were carried out to ensure the stability of these categories. The “proportion of ambiguous clustering” (PAC) was applied to infer the optimal number of clusters, where the K value has the lowest PAC.

The common genes in TCGA-PRAD DEGs and ERS-related genes were selected for univariate Cox regression analysis. Then, these BCR-related genes were retained for the LASSO model with the glmnet package, and 10-fold cross-validation was accepted to select the minimal penalty term. Then, the remaining genes were used to establish an optimal ER stress-related risk model, using the Akaike information criterion (AIC) method of multivariate Cox regression analysis. The ERS signature risk score was calculated as follows: risk_score = ∑ _{i =1} ⁿ (Coefi × Expi), where Coefi is the corresponding regression coefficient evaluated by multivariate Cox regression model and Expi is the expression value of the ERS-related genes. We divided patients into high-risk and low-risk groups based on the median risk score. Additionally, the TCGA-PRAD cohort was used as a training set to evaluate the prognostic value of BCR-dependent receiver operating characteristic (ROC) curves at 1, 3, and 5 years; GSE21034 and GSE70770 were utilized for the validation cohorts.

Drug-response prediction was evaluated based on the V2 database (809 cell lines and 198 compounds) of the Genomics of Drug Sensitivity in Cancer (GDSC) (15), using the “oncoPredict” package (16), and the half-maximum inhibitory concentration (IC₅₀) of each patient was assessed using the Ridge Regression model.

The R package “TCGAmutations” was used to calculate the total somatic mutation of TCGA-PRAD between different risk score subgroups. The online tool Sangerbox (17) (http://vip.sangerbox.com) was adopted to map an oncoplot to show the mutation landscape in high-risk and low-risk groups.

Based on the TCGA-PRAD dataset, a nomogram was established to predict the 1-, 3-, and 5-year BCR prognosis of PCa patients with ERS-related risk score and other related clinical parameters, using the “rms” package. Decision curve analysis and C-index were used to validate the clinical reliability of the nomogram model.

Univariate Cox and multivariate Cox regression were performed with risk score and other clinical variables to identify whether the ERS score was an independent prognostic predictor. The correlation of ERS-related genes’ risk score and age, Gleason score, and T stage were calculated.

The ESTIMATE package was used to evaluate the immune and stromal scores of PCa. The MCPcounter and CIBERSORT packages were used to detect the infiltration level of 22 immune cells and two stromal cells. Additionally, we also analyzed the tumor immune dysfunction and exclusion (TIDE) and MSI score for PCa patients in the two groups via TIDE tool (http://tide.dfci.harvard.edu). The immune checkpoints SIGLEC15, CTLA4, CD274, IDO1, PDCD1, HAVCR2, PDCD1LG2, and LAG3 between two groups were also analyzed.

The bioinformatics analysis was performed with R 4.1.3 (https://www.r-project.org/). Mann–Whitney–Wilcoxon tests or Student t-tests were carried out to analyze continuous variables. Survival plots were calculated by a log-rank test using the “survival” packages. The drug response was tested with Spearman’s correlation analysis. We chose p < 0.05 as the statistical significance.

To characterize the role of ERS-related genes, we extracted 787 ERS-related genes with correlation scores ≥7 from the GeneCards, as reported by Zhang et al. (12). We then analyzed the DEGs of these ERS-related genes in the TCGA cohort. As shown in Figure 1A , we obtain 108 DEGs. GO and KEGG enrichment analysis suggested that these DEGs were mainly enriched in the calcium ion homeostasis, calcium signaling pathway, glycosaminoglycan binding, PI3K-Akt signaling pathway, and ER lumen and adrenergic signaling in cardiomyocytes ( Figures 1B, C ). Additionally, based on the ERS-related DEGs, we classified PCa patients into two clusters with consensus clustering ( Figure 1D ). The cumulative density function and PAC method were used to identify the optimal k value ( Figures 1E, F ). We further performed survival analysis on these two clusters and found that the trend of BCR survival was worse in cluster 1 than in cluster 2 ( Figure 1G ), but not statistically significant, probably because of the number of queues, which requires further study.

To further construct the BCR prognostic model, we performed univariate Cox regression on ERS-related DEGs and identified 18 genes that were associated with BCR prognosis ( Figure 2A ). Furthermore, we analyzed these 18 genes by LASSO regression and derived 5 genes based on the minimum partial likelihood deviation ( Figures 2B, C ). Next, a five-gene BCR prognostic risk model was constructed by multivariate Cox regression analysis based on the AIC value ( Figure 2D ). Subsequently, the PCa cohort was divided into two subgroups based on the median risk score according to the risk model, and the expression of these five genes was shown in the heatmap ( Figure 2E ). BCR-free survival was analyzed using Kaplan–Meier and log-rank tests, and the results showed that the high-risk group had a shorter BCR-free time than the low-risk group ( Figure 2F ). In addition, the ROC curve was applied to assess the predictive efficiency of the model. The area under the ROC curve (AUC) was 0.775, 0.794, and 0.701 for 1, 3, and 5 years, respectively ( Figure 2G ), which indicated that the model performed well in predicting BCR-free survival in the TCGA cohort.

To better evaluate the predictive efficiency of the ERS-related gene model, we selected GEO datasets GSE70770 and GSE21034 for further validation. Risk scores were calculated for all patients in both datasets according to the same formula. We then divided the two datasets into high-risk and low-risk groups according to the median risk score. The edgeR package was then accepted to evaluate the expression patterns of these five genes in two GEO datasets and was found to be similar to TCGA ( Figures 3A, D ). As expected, survival time without BCR was significantly shorter in the high-risk group compared to the low-risk group ( Figures 3B, E ). Furthermore, we detected the predictive efficiency of this five-gene risk model utilizing ROC curves in GSE70770 and GSE21034, and the AUCs of 1, 3, and 5 years in both datasets suggested that the model was stable ( Figures 3C, F ).

To investigate the function of this ERS-related gene risk model in the clinical characteristics of PCa, we assessed the association between risk scores and clinical features. Results showed that compared with the low-risk group, the high-risk group had a higher tumor Gleason score and a more aggressive tumor stage ( Figures 4A, C ). The age of high-risk cohort patients was higher than that of the low-risk group ( Figure 4B ). In addition, univariate and multivariate Cox regression analysis suggested that the risk score was an independent prognostic factor in the TCGA dataset ( Figures 4D, E ) and extra GEO cohorts ( Figures S1D–G ). To better predict the BCR-free survival of PCa, we constructed a nomogram with ERS-related genes’ risk scores, age, tumor stage, and Gleason scores ( Figure 4F ). ROC curves were used to estimate the performance of nomogram, age, risk score, and tumor stage in predicting 1-, 3-, and 5-year BCR-free survival ( Figures S1A–C ). The calibration curves were utilized to assess the predictive performance of this nomogram with respect to the actual observed BCR-free survival rate ( Figure 4G ). The DCA curves indicated that the ERS-related gene risk signature had a favorable predictive efficiency ( Figure 4H ).

To identify the underlying mechanisms in different risk groups, we analyzed the DEGs in different risk subgroups in the TCGA cohort, and obtained 547 DEGs (|logFC| ≥ 1, p-value < 0.05) ( Figure 5A ). GO enrichment analysis was performed to assess the function of these DEGs, and the result showed that they were mainly involved in signaling receptor activator activity, receptor ligand activity, and hormone activity ( Figure 5B ). In addition, KEGG and GSEA showed that these ERS-related DEGs were mainly enriched in the IL-17 signaling pathway, ER protein processing, and TNF signaling pathway ( Figures 5C, D ). Furthermore, gene set variation analysis showed that the high-risk group was enriched in inflammatory and immune-related pathways, such as primary immunodeficiency and the Nod-like receptor signal pathway. Pathways related to histidine metabolism, phenylalanine metabolism, and calcium signaling were also varied in two groups ( Figure 5E ). The above results show that these DEGs participate in various processes including immune activity.

To further understand the potential characteristics of these two populations, we analyzed the top 15 individual cell mutation genes in both the high- and low-risk groups. TP53, TTN, SPOP, SYNE1, and KMT2D were the top five genes with the highest mutation frequency in the high-risk group, while in the low-risk group, they were SPOP, TTN, TP53, MUC16, and FOXA1 ( Figures 6A, B ). Meanwhile, patients in the high-risk group had higher levels of total mutational burden (TMB) than those in the low-risk group ( Figure 6C ). These results suggested that ERS-related genes may act through genetic mutations. Gene mutations can generate new antigens, and we wanted to know if the expression of immune checkpoints and immune infiltrating cells differed in the two tumor immune environments. We uncovered that PCa in the high-risk group had higher stromal and immune score ( Figures 6D, E ). MCPcounter showed that the high-risk group PCa had a higher abundance of T cells, NK cells, monocyte lines, and neutrophils ( Figure 6F ). CIBERSORT suggested that resting dendritic cells, Tregs, and macrophages M1 and M2 were significantly enriched in the low-risk group ( Figures 6G , S2 ). TIMER analysis showed that DNAJB1, COL10A1, PTGS2, AFP, and EGF were correlated with tumor-infiltrating lymphocytes ( Figure S3 ). Taken together, these data indicated that the different BCR-free prognosis may be related to the infiltrating immune cells.

Since the current study indicates that ERS-associated signature scores are associated with immune cell infiltration, we sought to find a link between ERS-associated gene risk scores and immune responses. This study evaluated the immune checkpoints in different subgroups based on risk score. The results showed that the expression levels of many immune checkpoints were lower in the low-risk subgroup than in the high-risk subgroup ( Figure 7A ). Moreover, the current study calculated TIDE scores between these two risk types, and the data showed that the low-risk group had a higher MSI score but a lower TIDE score and a lower T-cell dysfunction score ( Figures 7B–E ), indicating that the T-cell immune checkpoint inhibitor may be more effective in patients with a low-risk score. Furthermore, based on the GDSC database, the current study uses the “oncoPredict” package to look for compounds that may interact with ERS-related pathways. The results showed that ERS risk score-related drug sensitivity was related to apoptosis regulation (MIM1, WEHI-539, and ABT737), cell cycle (AZD7762), and WNT (WIKI4). ERS risk score-related drug resistance was associated with ERK MAPK (Selumetinib, SCH772984 and PD0325901), PI3K/MTOR (AZD2014), and EGFR (Gefitinib) ( Figures 7F, G ).