Transcriptional data (read counts), clinical characteristics and somatic mutation data were acquired from the TCGA database (Supplementary Table S1). Next, we downloaded four datasets with the same platform (Affymetrix GPL570) from the GEO database: GSE69223 (Meller et al., 2016) (N = 30), GSE32448 (Derosa et al., 2012) (N = 80), GSE55945 (Arredouani et al., 2009) (N = 21), and GSE46602 (Mortensen et al., 2015) (N = 50). Thereafter, we adjusted the background by using the “RMA” algorithm of the “affy” R package (Gautier et al., 2004) and removed the batch effect by the “ComBat” algorithm of the “sva” package (Leek et al., 2012). Therefore, we can merge the four datasets as the validation cohort. Moreover, GSE21034 (Taylor et al., 2010) (N = 370) was utilized to validate the prognostic value.

COSMIC is currently the broadest database of mutations in human cancer (Forbes et al., 2017). COSMIC mutation data (Genome Screens) and corresponding sample features were downloaded, and we estimated the frequency of each mutation site. Chi-square test was applied to discover the mutated genes with significant differences between primary and metastatic prostate cancer. Next, we executed prognostic analysis for each gene discovered by univariate Cox regression, and the genes related to prognosis with p-value < 0.05 were extracted for further analysis.

Unsupervised consensus clustering of the obtained genes was executed by using the k-means algorithm in the “ConsensusClusterPlus” package (Wilkerson and Hayes, 2010), which was repeated 1,000 times to ensure the stability of the classification. Survival differences between the two clusters were visualized by Kaplan–Meier curves. To explore the molecular characteristics of the two clusters. The “c2.cp.kegg.v7.4.symbols” gene set was downloaded from the MSigDB database, and we applied the “GSVA” package (Hänzelmann et al., 2013) to perform the GSVA analysis.

The “ssGSEA” method was performed to estimate the infiltration degrees of 28 immune cells by using the “GSVA” package. Estimate is commonly used to calculate scores reflecting the infiltration levels of immune cells and stromal cells in the tumor microenvironment by the package “estimate” (Yoshihara et al., 2013). We applied the estimate algorithm to calculate the ImmuneScore and StromalScore of each sample. Additionally, the correlation between the expression of immune checkpoint genes and androgen receptor between the two clusters was estimated.

The Pearson correlation coefficients of mutated genes with the identified PCa subclasses were estimated, and the signature genes A and B were obtained based on the correlation coefficients. Then, we applied the Boruta algorithm to the positively and negatively correlated genes to select feather genes. Finally, principal component analysis (PCA) was performed to estimate the first principal components of signature genes A and B. We defined the riskscore of each sample as follows: Riskscore = ∑ P C 1 B − ∑ P C 1 A

The difference in the riskscore in patients stratified by clinical parameters was estimated to expound the effect of the riskscore on cancer progression. Moreover, immune cell infiltration, ImmuneScore, StromalScore and the expression of immune checkpoint genes were also assessed between the high- and low-risk groups.

Tumor mutation burden (TMB) has been demonstrated as a predictive biomarker to identify whether patients with cancer can respond to immune checkpoint inhibitors well (Merino et al., 2020). Here, we explored the correlation between TMB and riskscore. Furthermore, we divided patients into four subgroups based on the median value of riskscore and TMB. Survival differences of the four subgroups were visualized by Kaplan–Meier curves.

Since the comparison of the expression of different immune checkpoint genes between the high- and low-risk groups was performed, here we used an immunotherapeutic cohort (IMvigor210 cohort) as a validation cohort (Mariathasan et al., 2018). We first evaluated the influence of the riskscore on the prognosis of patients treated with immunotherapy. Then, the riskscore of patients with different clinical statuses after treatment were compared. Finally, transcriptional data of tumor cell lines and IC50 values of antitumor drugs from the GDSC database were used to perform the drug sensitivity analysis by using the “pRRophetic” package (Geeleher et al., 2014).

All analyses were performed in RStudio 4.0.4. Correlation analysis was computed by the Spearman method. Student’s t test and the Wilcoxon test were applied for two-group comparisons. Correspondingly, the Kruskal–Wallis test and one-way ANOVA were used for multiple groups. Statistical significance was defined as p-value < 0.05.

The roadmap of this study is illustrated in Figure 1. The top 20 mutated genes for prostate cancer are shown in Figure 2A, and the mutation frequencies are displayed next to the gene name. Furthermore, Figures 2B,C shows the distribution of different types of mutations for prostate cancer. Missense substitution (88.07%), synonymous substitution (49.48%) and nonsense substitution (37.91%) were the main types of mutations, and the substitution mutations mainly included G > A (72.93%), C > T (72.42%), A > G (64.57%), and G > T (61.27%). Moreover, the Manhattan plot depicted mutation sites that had significant differences between primary and metastatic prostate cancer (Figure 2D).

We obtained 3,716 mutated genes with significant differences between primary and metastatic prostate cancer by the Chi-square test (Supplementary Table S2). Thereafter, we explored the prognostic value of these genes for progression-free survival (PFS) by the univariate Cox method, and 1,051 genes were extracted (Supplementary Table S3) for further analysis.

After comprehensive consideration of CDF curves and delta area, we chose k = 2 as the optimal cluster number for the clustering (Figure 3). Finally, 308 patients were assigned to Cluster 1, and 187 patients were assigned to Cluster 2. We also applied t-SNE dimension reduction, and the results suggested that the discrimination among subgroups was decent (Supplementary Figure S1). Survival curves suggested that Cluster 2 had a significant survival advantage compared with Cluster 1 (Figure 3C). Moreover, GSVA analysis and limma analysis (log FC > 0.2, adjusted p-value < 0.05) were performed. Significant differences in pathways related to cancer progression, such as the ERBB signaling pathway and VEGF signaling pathway, and pathways associated with the immune response, such as the B cell receptor signaling pathway and T cell receptor signaling pathway, were observed between the two clusters (Figure 3D).

In order to delve into the immune-related characteristics of the two clusters, the infiltration levels of immune cells were compared between the two clusters. A significant difference was observed in the infiltration degree of all immune cells, and all immune cells infiltration were lower in Cluster 2 (Figure 4A). Moreover, the results indicated that both the StromalScore and ImmuneScore of Cluster 2 were significantly lower compared with Cluster 1 (Figure 4B). What’s more, the expression of immune checkpoint genes, including CTLA4, PD-1, PD-L1 and PD-L2, appeared to be decreased in Cluster 2 (Figures 4C–F). However, the expression of AR was higher in Cluster 2 than in Cluster 1 (Figure 4G), which is consistent with the poor survival of Cluster 2.

Previous results indicated that the subclass was closely associated with the prognosis and immune infiltration levels of patients. However, this population-based classification cannot be directly used in clinical practice. Therefore, we constructed a scoring system to estimate the riskscore to predict the outcome of the patients and aid in making treatment decisions. After performing the Boruta algorithm, 230 genes that were positively and negatively correlated to the subclass were defined as signature genes A and B (Supplementary Table S4). The riskscore was acquired by conducting PCA on each signature gene (Supplementary Table S5). Patients were classified into high- and low-risk groups according to the cut-off point gained by the “survminer” package (Supplementary Figure S2). The results of survival analysis showed that patients in the high-risk group had lower PFS than those in the low-risk group (Figure 5A). As shown in Figure 5B, the high-risk group possessed a higher proportion of death. In the validation cohort GSE21034, patients with higher riskscore also showed a significantly shorter median PFS (Figure 5C). Moreover, it was observed that the riskscore was elevated in the high-risk clinical group with the progression of tumor (Figure 5D), and patients who achieved complete response after treatment had a significantly lower riskscore than other outcomes (Supplementary Figure S3).

TMB is emerging as a potential biomarker to predict the patients response to immune checkpoint inhibitors (Chan et al., 2019). Here, we evaluated the association between the riskscore and TMB. As shown in Figure 6A, patients in the high-risk group had a higher TMB than those in the low-risk group, and patients with a higher TMB had lower PFS (Figure 6B). Moreover, the correlation analysis indicated that the riskscore was positively associated with TMB (Figure 6C). Next, we combined the riskscore and TMB to divide patients into four subgroups. The patients with a high riskscore and high TMB had the shortest median PFS, and patients with a low riskscore and low TMB performed the best prognosis (Figure 6D). Finally, the mutation status of genes with high mutation frequencies in the high- and low-risk groups was visualized (Supplementary Figure S4).

Anticancer immunotherapies involving immune checkpoint inhibitors have emerged as new therapeutic regimens (O’Donnell et al., 2019). The tumor microenvironment (TME) was proven to be tightly linked to tumor progression and metastasis (Brassart-Pasco et al., 2020) and can blunt the therapeutic response, thus affecting the clinical outcome (Wu and Dai, 2017). To further explore the correlations between patients’ response to immunotherapy and riskscore, we first compared the immune cell infiltration levels between high- and low-risk groups, and the results indicated that compared with the low-risk patients, the infiltration levels of 28 immune cells in the high-risk patients were significantly downregulated (Figure 7A).

Moreover, infiltration estimation for TCGA was downloaded from TIMER 2.0 (Li et al., 2020), which includes several methods, such as XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, and CIBERSORT. We correlated the immune cell infiltration levels and the riskscore, and the results also suggested that the riskscore was negatively related to most of the immune cells (Supplementary Figure S5). Additionally, the StromalScore and ImmuneScore of the high-risk patients were significantly lower than low-risk patients as well (Figure 7B). Furthermore, the associations between immune checkpoint inhibitor genes and the riskscore were evaluated. The expression of ICI genes, such as CTLA4, PD-1, PD-L1 and LAG3, were downregulated in the high-risk groups compared with the low-risk groups (Figure 7C). Finally, we validated the performance of the riskscore in the IMvigor210 cohort. As shown in Figure 8A, the low-risk patients still showed a significant survival advantage compared with high-risk patients. We evaluated the differences in riskscore among patients with different responses to immunotherapy. Patients with complete response (CR) had the lowest riskscore while patients performed progressive disease (PD) had the highest riskscore (Figure 8B). These results suggested that patients in the high-risk group may not have a satisfying response to immunotherapy.

Considering that androgen deprivation therapy and chemotherapy still play vital roles in the treatment of prostate cancer. We used the GDSC database to explore the association between the riskscore and drug sensitivity. The results revealed little difference in the predicted IC50 of Bicalutamide between the high- and low-risk groups (Figures 8C,E). However, both in the TCGA and GEO cohorts, significant differences were observed in the predicted IC50 of Docetaxel between the high- and low-risk groups (Figures 8D,F). The IC50 of Docetaxel was significantly lower in the high-risk group, suggesting that these patients are sensitive to Docetaxel.