Raw counts of RNA-seq data and clinical data of 144 PRAD patients were obtained from International Cancer Genome Consortium (ICGC https://www.icgc-argo.org), and another 499 samples from The Cancer Genome Atlas (TCGA https://www.cancer.gov/tcga). Then the RNA-seq counts were transformed into transcript per million (TPM) to avoid errors in data processing. A total of 2,483 immune-related symbols were retrieved from the Immport Database (Bhattacharya et al., 2018) (https://browser.immport.org/browser). Finally, 2026 immune-related genes were identified for follow-up data analysis after comparison with Req-seq data from the TCGA database.

Gene Expression Profiling Interactive Analysis (Tang et al., 2019) (GEPIA2, http://gepia2.cancer-pku.cn, version 2) was used to calculate the differential gene expression and patient survival data. One-way ANOVA was used to identify differentially expressed genes with q values < 0.05 and |log2FC| > 1. The Kaplan-Meier method was used to evaluate OS and PFS, with a cut-off value of 50% (median), and the log-rank test was used for comparison. And adjusted p-value < 0.05 was considered statistically significant.

We used cBio Cancer Genomics Portal (Cerami et al., 2012) (cBioPortal, http://www.cbioportal.org) to analyze copy number alterations (CNAs), mutations in the raw RNA-seq data from TCGA databases and compare genetic alterations in PRAD.

The relationship between the abundance of tumor immune infiltrating cells (TIIC) and mRNA vaccine antigens was analyzed on Tumor Immune Estimation Resource (Li et al., 2017) (TIMER, https://cistrome.shinyapps.io/timer/). Spearman correlation analysis was used to select the purity adjustment. p-value < 0.05 was considered statistically significant.

Clusters were made based on their gene expression profiles of immune cell lineages. The 2026 immune-related genes were clustered and a consistency matrix was established to identify corresponding gene modules and immune subtypes. The partition around medoids (PAM) algorithm was applied with distance measured as ‘1—Pearson’ correlation coefficient, and 500 replicates of bootstraps were performed, each involving 80% patients of the TCGA cohort. K value of cluster varied from 2 to 6 and the optimal partition was determined by evaluating the consensus matrix and the consensus cumulative distribution function. We got consistent results in the ICGC cohort in the same settings. The consistency of immune subtypes between TCGA and ICGC cohorts was quantified by calculating the Pearson correlation between the proportion within the group and the centroid of gene module score.

The log-rank test was used to evaluate the prognostic values of the immune subtypes, and univariate Cox regression analysis was performed with PFS as the endpoint. The association of immune subtypes with different immune-related molecular and cellular characteristics was determined using the Kruskal–Wallis test.

To visualize the distribution of immune subtypes across individual patients, the reduceDimension function was used to perform dimension reduction with Monocle implementation of Discriminative Dimension Reduction Tree (DDRTree) algorithm (Qi et al., 2017). The maximum number of components was set to 5. Finally, the immune landscapes were displayed with color-coded immune subtypes using the function plot cell trajectory.

The genes intersecting between LASSO and Boruta analysis were considered the most critical relevant variables for further analysis, which were visualized by a Venn diagram. Then, the discriminative performance of the predictors was evaluated by receiver operating characteristic (ROC) curve analysis, and the optimal cutoff values, AUCs, sensitivity, specificity, and accuracy were determined.

To identify potential antigens of PRAD, 18,335 amplified genes (Supplementary Table S1) and 3,016 overexpressed genes (Supplementary Table S2) that could express the tumor-associated antigens were explored (Figures 1A,B). Then, a total of 10,402 mutant genes (Supplementary Table S3) potentially encode tumor-specific antigens were screened by evaluating fraction genome alteration and mutation counts in individual samples (Figures 1C,D). Genes with the highest mutation frequency in the mutation count group, including ETS transcription factor ERG, Transmembrane serine protease 2, Phosphatase and Tensin homolog, Titin, STAM binding protein-like 1, Lipase family member N, Lipase family member M, Lipase family member K, Lipase family member F and Ankyrin repeat domain 22, were individually displayed (Figure 1E). Gene with the highest mutation frequency in the fraction genome altered group, including LDL receptor-related protein 1B, ETS transcription factor ERG, Transmembrane serine protease 2, ST3 beta-galactoside alpha-2,3-sialyltransferase 1, thrombospondin type 1 domain containing 7B, Zinc finger protein 292, ring finger protein 19A, RBR E3 ubiquitin-protein ligase, vacuolar protein sorting 13 homolog B and MYC proto-oncogene, were displayed separately (Figure 1F). Consistent with previous reports, the ERG, TMPRSS2, PTEN, and zinc finger protein family are closely related to the occurrence and progression of prostate adenocarcinoma. For example, TMPRSS2-ERG gene fusion often predicts the progression and poor prognosis of prostate adenocarcinoma (Linn et al., 2016). Phosphatase and tensin homolog (PTEN) loss has long been associated with adverse findings in early prostate adenocarcinoma (Jamaspishvili et al., 2020). In total, 1,003 genes with amplification, overexpression, and mutation were identified.

A total of 40 genes associated with the OS of PRAD were identified, of which 7 genes were related to the PFS (Figure 2A). As shown in Figure 2B, the univariate Cox regression analysis showed that the expression of the above 7 genes was related to the prognosis of PRAD patients. Figures 2C,D showed that the chromosomal distribution and genomic variation of 7 genes above-mentioned. The elevated expression of ADA (Figures 2E,F), FYN (Figures 2G,H), HDC (Figures 2I,J), NFKBIZ (Figures 2K,L), RASSF4 (Figures 2M,N), SLC16A3 (Figures 2O,P) and UPP1 (Figures 2Q,R) indicated a worse prognosis of PRAD. Furthermore, the higher expression levels of ADA (Figure 3A), FYN (Figure 3B), HDC (Figure 3C), NFKBIZ (Figure 3D), RASSF4 (Figure 3E), SLC16A3 (Figure 3F), and UPP1 (Figure 3G) were significantly correlated with increased tumor infiltration of macrophages, Dendritic cells (DCs), and B cells, suggesting that the seven tumor antigens may have potential immune-stimulatory effects. In summary, a total of 7 candidate genes with critical roles in PRAD development and progression were identified. These findings indicate that seven tumor antigens (ADA, FYN, HDC, NFKBIZ, RASSF4, SLC6A3, and UPP1) have potential immune-stimulating effects and can be processed and presented by antigen-presenting cells (APC) to induce tumor responses, and are hopeful candidates for developing anti-PRAD mRNA vaccine.

Through immunotyping, we can recognize the immune status of tumor and its microenvironment, and determine the patients suitable for vaccination. Then we analyzed the expression data of 2026 immune-related genes in 499 PRAD samples from the TCGA cohort to construct consensus clustering. Based on their cumulative distribution function and function δ area, we chosen k = 5, where immune-related genes seem to be stably clustered (Figures 4A,B), and 5 immune subtypes designated as C1-C5 were obtained (Figure 4C). C1 and C2 showed a significant correlation with superior prognosis, whereas C3 showed the poorest survival probability (Figure 4D). The distribution of different T stages and Gleason score in subtypes showed that patients diagnosed in different stages were clustered dispersedly. The Gleason score 10 group and stage T4 were significantly related to C3 (Figures 4E,F). Consistent with the results obtained from the TCGA cohort, the Gleason score 10 group also showed a significant correlation with C3 in the ICGC cohort (Figure 4G). In conclusion, immunotyping can be used to predict the prognosis of patients with PRAD, and its accuracy is better than traditional Gleason score and staging, which is consistent across different cohorts.

Higher tumor mutational burden (TMB) and somatic mutation rates are associated with stronger anticancer immunity (Rooney et al., 2015). Therefore, we calculated the TMB value and mutation for each patient and performed the same analysis in all immune subtypes using the TCGA mutation data set processed by mutect2. C1 and C2 showed much higher TMB and mutated genes counts compared to C5 (Figures 5A,B). Furthermore, Comparing with C5 (53.01%), C1 and C2 had higher mutation rates (81.31%, 81.82%). Surprisingly, C3 and C4 also showed higher mutation rates (81.82%, 63.46%) among the five subtypes (Figures 5C–G). The mutation frequency of 11 famous PRAD related genes including SPOP, TTN, TP53, MUC16, and FOXA1, etc. were the highest in each subtype (Figure 5H). In all, our findings indicate that immunotyping can be used to predict TMB and somatic mutation rates in PRAD patients and that patients in type C1 and C2 with higher TMB and mutated genes counts may have positive feedback on mRNA vaccines.

Immune checkpoints (ICPs) and immunogenic cell death (ICD) modulators play critical roles in regulating host anti-tumor immunity. We then analyzed their expression levels in different immune subtypes. A total of 47 ICPs related genes were analyzed in TCGA and ICGC cohorts, of which 46 (97.9%) genes in the TCGA cohort (Figure 6A) and 43 (91.5%) in the ICGC cohort (Figure 6B) had differential expression between the immune subtypes. For example, CD200, CD200R1, CD244, CD27, CD28, CD40, CD40LG, CD48, CD80, CD86, CTLA4, HAVCR2 and ICOS, etc. were overexpressed in the C3 and C5 tumors in the TCGA cohort, which is consistent in the ICGC cohort.

In addition, the general expression level of ICPs in the TCGA cohort was higher than that in the ICGC cohort. 26 ICD related genes were detected in the TCGA cohort, of which 23 (88.5%) had differential expression between the immune subtypes (Figure 6C). Similarly, 26 ICD related genes were detected in the ICGC cohort, of which 21 (80.8%) had differential expression between the immune subtypes (Figure 6D). For example, CALR, IFNAR1, EIF2AK3, EIF2AK2, EIF2AK4, and EIF2AK1 showed higher expression levels in C1 and C2 tumors in both cohorts. Collectively, immunotyping can mirror the expression level of ICD modulators and ICPs, and be used as potential therapeutic biomarkers for mRNA vaccines. According to the above results, C3 and C5 were more likely to show immunosuppressive phenotype. By comparison, patients with C1 and C2 are suitable for mRNA vaccine.

Gleason score and serum PSA level are currently recognized indicators for the prognosis and diagnosis of PRAD. The higher levels of them often indicate the progress, poor prognosis, or recurrence of PRAD. In this study, the levels of Gleason score and serum PSA were analyzed for each patient in TCGA and ICGC cohorts. Both TCGA and ICGC cohorts showed significant differences in Gleason score and serum PSA level among the immune subtypes. For example, compared with C3 and C5, C1 and C2 had significantly lower Gleason score and serum PSA level in the TCGA cohort (Figures 7A,B). These results were consistent with the superior prognosis showed in the C1 and C2 patients. It is also confirmed that Gleason score and serum PSA level are important prognostic markers of PRAD. Taken together, in predicting the prognosis of PRAD patients, the immune subtype is not inferior to the Gleason score and serum PSA level.

The effectiveness of mRNA vaccines largely depends on the immune status of the tumor. Therefore, we used ssGSEA to score 28 previously reported characteristic genes in the TCGA and ICGC cohorts to define immune cell components in the five immune subtypes. As shown in Figure 8A, the immune cell components were divided into 5 subtypes. There were significant differences in immune cell composition among different subtypes. Almost all the immune cell related genes mentioned in C3 and C5 groups were significantly higher than those in C1 and C2 groups in TCGA cohorts (Figures 8A,B). Thus, C3 and C5 were immune “hot” and immunosuppressive phenotypes, while C1 and C2 were immune “cold” phenotypes. A similar trend occurred in the ICGC cohort, as well (Figures 8C,D). Thus, these results indicate that the immune subtype can mirror the PRAD immune status, and identify patients suitable for mRNA vaccination. The induction of immune responses in patients with immune “cold” C1 and C2 tumors treated with the mRNA vaccine seems promising.

To prove the reliability of this immune classification, we next explored the correlation between the five immune subtypes and the six previously reported immune classifications (Thorsson et al., 2018). For the convenience of clear description, they were named IS1-IS6. Of which PRAD was mostly clustered into IS3. C1 and C2 mainly overlapped with IS3 and IS4, C3 with IS3 and IS1, and C5 with IS1, IS2 and IS3 (Figure 8E). IS3 was associated with superior and worse prognoses respectively, IS1 and IS2 indicated intermediate prognoses. These findings were consistent with the observed longer survival of patients with C1and C2 tumors than those with C3 and C5 tumors. Surprisingly, C1 patients with superior prognosis and C3 patients with the poorest survival largely overlapped with IS3. These results not only proved the reliability of our immunotyping method but also augmented the previous classification. We also tried to obtain 28 immune cell markers from the public database and used the ssGSEA algorithm to calculate the immune cell infiltration rate in TCGA and ICGC cohorts. The correlation between the immune subtypes and identified 28 significantly associated immune cells, as well as 17 associated immune signatures in the TCGA cohort, were analyzed (Charoentong et al., 2017). As shown in Figures 8F,G, C3 and C5 had the most obvious immune characteristics and much higher immune cell infiltration rates in almost all groups. Therefore, C3 and C5 were associated with an overall favorable immune-activated phenotype and characterized by diverse immune signatures. In contrast, the lower scores of immune signatures in C1 and C2 were indicative of an immune ‘cold’ phenotype. To sum up, immune subtypes can reflect the cellular and molecular characteristics of patients with PRAD and also indicate their immune status. Therefore, immune subtypes are hopeful biomarkers for mRNA vaccines. Patients without the immunosuppressive microenvironment, such as immune “cold” C1 and C2 tumors, may be suitable candidates for mRNA vaccination.

In this study, we constructed the immune landscape of PRAD using the immune gene expression profiles of individual patients (Figure 9A). The X-axis correlated with a variety of immune cells, of which plasmacytoid dendritic cells and CD56^dim natural killer cells showed the most negative correlation, and the Y-axis was most negatively correlated with CD56^dim natural killer cells (Figure 9B). The integral distribution of C3 and C5 was opposite to that of C1 and C2. Surprisingly, even the same subtypes show the opposite distribution, which indicates significant intra-cluster heterogeneity within the subtypes, especially in C4 and C5. C4 was further divided into three subsets and C5 into two subsets according to the location of the immune cell population (Figure 9C). There were significant differences in the enrichment scores of several immune cells between the subtypes. (Figure 9D). For instance, C4A and C5A had lower counts of activated B cells, activated CD4 + T cells, activated CD8 + T cells, effector memory CD8 + T cells, regulatory T cells, and bone marrow-derived suppressor cells (MDSC). Therefore, mRNA vaccines may be relatively feasible and more effective in C4A and C5A. In addition, after comparing the prognosis of samples with extreme distribution in the immune landscape, we found that group 1 showed the best survival probability. This result is consistent with the above results (Figures 9E,F). In conclusion, these findings suggest that the immune landscape based on immune subtypes can accurately identify the immune components of each PRAD patient and predict their prognoses, which is helpful to select the appropriate patients for mRNA vaccine.

Then, differential gene expression analyses were performed between C1 and C2 and other groups (Figure 10A). The most representative genes of C1 and C2 were used to identify immune “cold” tumor-related features using machine learning methods. A total of 12 genes (Supplementary Table S4) and 62 genes (Supplementary Table S5) were identified by least absolute shrinkage and selection operator (LASSO), binomial deviance, and Boruta analyses (Figures 10B–D) respectively, and the 12 genes shared by both methods were determined as the specific features of C1 and C2 tumors (Figure 10F). In addition, the C1 and C2 predictor genes also had an AUC of 0.919 and 0.828 in the TCGA cohort and the ICGC cohort (Figures 10G,H), respectively. Taken together, the gene-set of 12 predictors (IGF2, H2AFJ, DCXR, MCRIP2, CFAP65, ANTKMT, HPN, SCAND1, SLC25A10, RHPN1, MIF, and ZNF593) has an excellent performance in the identification of patients suitable for mRNA vaccine.