Inference of immune infiltration cells in TME
The landscape of TME cell infiltration models was inferred by the CIBERSOFT algorithm. Figure 1A,B summarizes the findings obtained from the 64 normal samples, 32 adjacent tumor samples, and 351 tumor samples. Twenty-one types of TIICs were detected in patients, while naïve CD4 + T cells did not appear in any samples, the findings were in accordance with the previous study[25]. The Correlation matrix of all 21 immune cell densities in the TCGA cohort was also shown in Fig. 1C,D. The results in current works found that the proportion of TME cells in PCa was significantly different among the three groups. In general, the infiltration of TME cells differed among these samples, with the highest proportion of LM22 immune cells found in normal samples, followed by adjacent tumor samples, and finally tumor samples. The mutual interaction relationship between immune cells in normal and adjacent tumor samples is evident compared to tumor samples (Fig. 2A,B,C,D). Compare to normal tissues, the fraction of total T cells and total B cells were higher in adjacent tumor tissues and tumor tissues. Total Macrophage was mainly observed in the tumor tissue (Fig. 3A,B,C). B cell memory was devoid in the adjacent tumor samples and tumor samples, Macrophages M1 and M2 were increased in the tumor samples, T cells gamma delta and T cells regulatory (Tregs) was found in the adjacent samples and tumor samples. We can infer from the results that the immune system is inhibited during tumor development.
We found that B cell memory, Macrophages M2, Mast cells activated, Mast cells resting, Monocytes, NK cells resting, NK cells activated, T cells CD8 and T cell follicular helper were mainly enriched in the normal tissues, while B cell naive and Dendritic cells activated, Dendritic cells resting, Macrophages M1, T cells CD4 memory activated, T cells CD4 memory resting, Neutrophils were significantly founded in the adjacent tissues. In the tumor tissues, B cell naïve, Dendritic cells activated, Dendritic cells resting, Macrophages M0, Macrophages M1, T cells CD4 memory activated, T cells CD4 memory resting and T cells regulatory (Tregs) were mainly activated (Fig. 4). Therefore, these results demonstrated that aberrant immune infiltration and it's heterogeneous in PCa as a tightly regulated process might play important roles in the development of tumor, it has important clinical meanings.
Prognostic Subsets Of Immune Cells
The univariate Cox regression analysis was performed to identify the prognostic subtypes of TIICs in PCa, and the results were shown that B cell memory was significantly correlated with RFS with the cut-off of the P-value less than 0.05 (Table 1, Fig. 5A). In the multivariate Cox regression analysis, B cell memory was still obviously related to RFS (Table 2, Fig. 5B). Then, we performed a Kaplan-Meier curve plot and log-rank test for the above-mentioned immune cell subsets, and the results are shown in Fig. 5C,D. B cell memory and T cells regulatory (Tregs) were significantly associated with RFS with PCa patients.
Table 1
The univariate Cox regression analysis was performed to identify the prognostic subtypes of TIICs in PCa.
Immune cell | HR | HR.95L | HR.95H | pvalue |
B cells naive | 0.663123 | 0.000441 | 996.5155 | 0.912357 |
B cells memory | 2.98E + 34 | 2.19E + 13 | 4.06E + 55 | 0.001388 |
Plasma cells | 0.001334 | 1.20E-08 | 148.7424 | 0.264269 |
T cells CD8 | 2.075994 | 0.005774 | 746.3596 | 0.80779 |
T cells CD4 memory resting | 0.253485 | 0.000796 | 80.75887 | 0.640723 |
T cells CD4 memory activated | 0.003227 | 2.03E-12 | 5119968 | 0.595629 |
T cells follicular helper | 0.000369 | 7.62E-10 | 178.984 | 0.236668 |
T cells regulatory (Tregs) | 79.00902 | 2.48E-08 | 2.52E + 11 | 0.695531 |
T cells gamma delta | 0.000811 | 3.25E-22 | 2.02E + 15 | 0.741926 |
NK cells resting | 2831.234 | 1.18E-30 | 6.77E + 36 | 0.839372 |
NK cells activated | 2.001437 | 2.76E-07 | 14491266 | 0.931388 |
Monocytes | 3357.781 | 0.000786 | 1.43E + 10 | 0.297283 |
Macrophages M0 | 11.70166 | 0.01414 | 9683.931 | 0.473023 |
Macrophages M1 | 2762.816 | 0.007955 | 9.6E + 08 | 0.223474 |
Macrophages M2 | 70.26358 | 0.016501 | 299185.2 | 0.318604 |
Dendritic cells resting | 0.056424 | 3.43E-06 | 927.6329 | 0.561617 |
Dendritic cells activated | 1.11E-14 | 3.10E-31 | 399.1275 | 0.098531 |
Mast cells resting | 2.98651 | 0.001924 | 4635.156 | 0.770392 |
Mast cells activated | 4502.05 | 1.40E-19 | 1.45E + 26 | 0.75037 |
Eosinophils | 8.19E-34 | 5.81E-99 | 1.16E + 32 | 0.319553 |
Neutrophils | 2.02E-39 | 3.99E-89 | 1.02E + 11 | 0.127062 |
Table 2
The multivariate Cox regression analysis was conducted to screen the prognostic subtypes of TIICs in PCa.
Immune cell | HR | HR.95L | HR.95H | pvalue |
B cells naive | 6.67E + 39 | 7.90E-12 | 5.63E + 90 | 0.125353 |
B cells memory | 6.60E + 69 | 9.74E + 12 | 4.48E + 126 | 0.016044 |
Plasma cells | 2.19E + 35 | 1.15E-16 | 4.18E + 86 | 0.176776 |
T cells CD8 | 3.96E + 37 | 7.90E-14 | 1.99E + 88 | 0.146099 |
T cells CD4 memory resting | 6.14E + 36 | 4.59E-15 | 8.20E + 87 | 0.158454 |
T cells CD4 memory activated | 2.66E + 37 | 5.70E-16 | 1.24E + 90 | 0.163714 |
T cells follicular helper | 1.92E + 34 | 3.94E-18 | 9.33E + 85 | 0.193608 |
T cells regulatory (Tregs) | 2.48E + 38 | 3.96E-14 | 1.55E + 90 | 0.146276 |
T cells gamma delta | 8.99E + 38 | 8.48E-17 | 9.52E + 93 | 0.165286 |
NK cells resting | 4.12E + 53 | 3.86E-18 | 4.40E + 124 | 0.139016 |
NK cells activated | 2.11E + 38 | 4.75E-13 | 9.35E + 88 | 0.138054 |
Monocytes | 2.04E + 44 | 3.07E-08 | 1.36E + 96 | 0.093777 |
Macrophages M0 | 1.19E + 38 | 2.38E-13 | 5.92E + 88 | 0.141039 |
Macrophages M1 | 4.64E + 40 | 1.24E-10 | 1.74E + 91 | 0.115014 |
Macrophages M2 | 2.82E + 38 | 3.46E-13 | 2.29E + 89 | 0.138805 |
Dendritic cells resting | 1.98E + 38 | 1.84E-13 | 2.13E + 89 | 0.141324 |
Dendritic cells activated | 3.37E + 20 | 4.47E-34 | 2.54E + 74 | 0.455198 |
Mast cells resting | 3.10E + 38 | 2.19E-13 | 4.37E + 89 | 0.140241 |
Mast cells activated | 2.80E + 55 | 0.00022 | 3.57E + 114 | 0.065967 |
Eosinophils | 4.63E-10 | 4.73E-95 | 4.52E + 75 | 0.82956 |
Neutrophils | NA | NA | NA | NA |
Patterns Of Immune Cell Infiltration In Molecular PCa Subclasses
The molecular classification of PCa was identified by unsupervised consensus clustering in all tumor samples using“CancerSubtypes” R-package. The optimal number of clusters is determined by the K value. After evaluating the relative changes in the area under the cumulative distribution function (CDF) curve and the consensus heat map, a three-cluster solution (K = 3) with no significant increase in area under the CDF curve was selected (Fig. 6). This result revealed 88 patients (38%) in Cluster I, 92 patients (48%) in Cluster II, and 133 patients (64%) in Cluster III for the PCa cohort. The consensus matrix heat map demonstrated that cluster I and cluster II, and clusters III appear well in the individualized clusters in Fig. 7(A, B, C). Patients classified as Cluster III had a good prognosis compared to Cluster I and Cluster II (P < 0.001, log-rank test).
Differentially expressed genes and immune cell patterns of Cluster I, Cluster II and Cluster III subtypes
The molecular subtypes of PCa cause differences in clusters I, II, and III subtypes, which exhibit activation of specific signaling pathways and different prognosis. The Kruskal-Wallis test was performed to identify subtype-specific genes/LM22 immune cells in each subtype. In terms of subtype-specific immune cells in PCa, Cluster I and III is defined by high levels of Dendritic cells resting compared with Cluster II. Cluster II was enriched by Macrophages M0 and NK cells activated; While Cluster III was defined by the high level of Dendritic cells activated, Dendritic cells resting, T cells CD4 memory activated and T cells regulatory (Tregs) (Fig. 8).
To uncover molecular differences among PCa molecular subtypes and identify subtype-specific biomarkers, unpaired Student's t-test was used to identify differentially expressed genes that were significantly associated with each subtype. Gene differentially expressed analysis was conducted with the cut-off point of FDR < 0.05. In subgroup I, a total of 294 immune related mRNAs (227 up-regulated genes and 67 down-regulated genes) were differentially expressed genes (DEGs) compared to other groups; In subgroup II, a total of 292 DEGs (76 up-regulated genes and 216 down-regulated genes) compared other groups; in subgroup III, a total of 162 DEGs (162 up-regulated genes and 95 down-regulated genes) were observed compared to other groups (Fig. 9).
In terms of clinical features among three clusters, Cluster I has a Higher Gleason score and PSA level compared to Cluster II and III. However, the PSA value and Age has no difference between the three Clusters (Fig. 10A,B,C). The heatmap illustrates the association of different clinical characters among the three subgroups (Fig. 10D). Statistical significance was performed by the Kruskal-Wallis test.
Identification of differentially expressed genes and enriched Gene Ontology and pathway in subtypes
Functional enrichment analysis for DEGs in CI vs CII, III, CII vs CI, III and C III vs CI, II was performed. For CI vs CII, III,a total of 517 GO terms of biological process, 25 GO terms of cellular component, and 39 GO terms of molecular function were identified with the cutoff point of adjust P-value < 0.05 (Supplementary Table 1). In subgroup II compared to other groups, 476 GO terms of biological process, 25 GO terms of cellular component, and 45 GO terms of molecular function were identified. In subgroup III compared to other groups, a total of 425 GO terms of biological process, 13 GO terms of cellular component, and 18 GO terms of molecular function were identified. Top GO terms included cytokine activity, immune/inflammatory response, and chemokine activity (Fig. 11A). All Clusters were enriched in cytokine activity, chemokine activity and carbohydrate binding. Cluster I and II were associated with lipopeptide binding, while Cluster III was related to receptor ligand activity and G protein − coupled purinergic nucleotide receptor activity, G protein − coupled nucleotide receptor activity (Fig. 12A).
Furthermore, all pathways generated by the Kyoto Genomics and Genomics Encyclopedia (KEGG) analysis are related to immune responses (Fig. 11B). Cluster I was associated with the Primary immunodeficiency, Th1 and Th2 cell differentiation and NF − kappa B signaling pathway, Cluster II was enriched for Primary immunodeficiency, JAK − STAT signaling pathway and Natural killer cell mediated cytotoxicity, while Cluster III was related to Th17 cell differentiation, Cell adhesion molecules (CAMs) and TNF signaling pathway (Fig. 12B).
To identify gene sets enriched in each subtype, we then performed GSEA analysis. GSEA analysis reveals distinct enriched gene sets between subtypes. The number of enriched pathways progressively increased from Cluster I through Cluster II. We subsequently selected representative gene sets for CI-CIII to build a pathway heatmap, which reveals distinct gene sets enriched in each subtype. Cluster I was associated with SMID BREAST CANCER ERBB2 and SHEDDEN LUNG CANCER POOR SURVIVAL, Cluster II was enriched for RUTELLA_RESPONSE_TO_HGF_VS_CSF2RB_AND_IL4_DN and VILIMAS_NOTCH1_TARGETS, while Cluster III was related to BASSO CD40 SIGNALING and VERHAAK AML WITH NPM1 MUTATED (Fig. 13A,B,C,D). All these similarities are in good agreement with the molecular and clinical characteristics of the three subtypes identified in our study, confirming the correctness of the features we found on these three subtypes.
Subclass-associated Gene Mutations
This study also investigated the association between the three subtypes and count of somatic mutation. Gene mutation profiles of these highly mutated genes are shown in Fig. 14A,B,C,D,E,F. SPOP showed a higher mutation rate in Cluster I, TP53 exhibited a higher mutation rate in Cluster II, while in the Cluster III, CACNA1E and TTN, TP53 exhibited a higher mutation rate.
Sensitivity to Immuno/Chemotherapies among Cluster I, II and III Subtype
Although immunological checkpoint drugs have not been approved for use as a conventional drug for PCa, we have used the TIDE algorithm to predict the likelihood of responding to immunotherapy, and it has demonstrated that cluster III was more likely better to the response to immunotherapy than clustering II and III (P < 0.05). As for the TIDE prediction, a subclass mapping method was employed to compare the expression profiles of the three PCa subtypes with another published data set containing 47 melanoma patients who responded to immunotherapy. The results showed that Cluster I was more promising for CTLA4 treatment (Bonferroni correction P < 0.05) (Fig. 15A).
Given that chemotherapy is currently a common treatment for prostate cancer, we sought to evaluate the response of three subtypes to common chemotherapy drugs. Therefore, we trained the prediction model on the GDSC cell line dataset by ridge regression and evaluated the satisfactory prediction accuracy by 10-fold cross-validation. We estimated the IC 50 for each sample in the TCGA dataset based on the predictive model of these two chemo drugs. For several common chemicals, we can observe that Cluster I and Cluster II, Cluster III have significant differences in estimated IC 50, of which Cluster I may be more sensitive to commonly used chemotherapy (AKT.inhibitor.VIII P < 0.001, ATRA P = < 0.001) (Fig. 15B).