Normalized RNA sequencing data (fragments per kilobase million, FPKM), based on mRNA samples (498 prostate cancer tissues and 52 normal tissues) from the Illumina HiSeq RNA-Seq platform of the Cancer Genome Atlas Prostatic Adenocarcinoma (TCGA-PRAD) database, were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) (Weinstein et al., 2013). FPKM data were transformed into transcript per million (TPM) values following log₂ (x + 1) normalization. The corresponding clinical profiles (age, AJCC-TNM stage, and Gleason score) were downloaded from the TCGA portal, and the clinical data for progression-free survival (PFS) analysis was downloaded from TCGA Pan-Cancer Clinical Data Resource (TCGA-CDR) (Liu et al., 2018). After removing four patients without survival data, 494 patients with PRAD were selected for the present study. The clinical features of the patients are summarized in Table 1.

The microarray data of 280 PRAD samples in the GSE116918 dataset, based on GPL25318 ([ADXPCv1a520642] Almac Diagnostics Prostate Disease Specific Array) (Affymetrix/Thermo Fisher, Belfast, United Kingdom), were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Detailed clinicopathological data included age, AJCC-T stage, Gleason score, metastasis-free survival (MFS) state, and MFS time. The expression data in the two platforms underwent a batch calibration process and were further normalized via the “sva” R package, such that they were comparable. The clinical features of the patients are summarized in Table 1. In addition, somatic mutation data of 475 patients with PRAD (based on the VarScan software) were obtained from the TCGA portal.

The 22 TIICs constituting the TME of TCGA-PRAD samples were calculated using the CIBERSORT algorithm (http://cibersort.stanford.edu/) (Newman et al., 2015). Samples with a CIBERSORT p-value <0.05 were used for further study.

The WGCNA R package was used to perform WGCNA of 19,560 gene sequences from TCGA-PRAD patients (Langfelder and Horvath, 2008). SampleTree was used to identify the outliers that were subsequently deleted. According to the mean connectivity and scale-free topology model fit, the soft threshold power (β) value was selected to generate a scaleless network (index of scale-free topologies = 0.90). The correlations between sample traits and candidate modules were computed to determine the models that were highly correlated with the traits. Then, similar genes were introduced into the same candidate module employing a “dynamic tree cutting” algorithm with a minimum size of 60. Correlations between the 22 TIICs and module characteristic genes were evaluated using Pearson’s correlation coefficient (p < 0.05). Finally, genes in the most statistically significant module were selected for subsequent analysis.

The expression levels of genes in the most statistically significant module were extracted from TCGA-PRAD and GSE116918 datasets. TCGA-PRAD was used as the training cohort, and GSE116918 was used as the validation cohort. Univariate Cox regression analysis was applied to obtain prognostic risk candidate genes from the most significant module in the “activated dendritic cells” population in the training cohort, and the genes that were significantly related to progression-free survival (PFS) (p < 0.01) were identified.

Then, least absolute shrinkage and selection operator (LASSO) Cox regression analysis was used to determine the best weighting coefficient of the prognostic risk candidate genes. After a 1,000-fold cross-validation of the maximum likelihood estimate of the penalty, the minimum criterion was determined using the optimal value of the penalty parameter λ. Finally, a TIIC-related gene risk signature was established, and the risk score was calculated using the following formula: R i s k   s c o r e = β   g e n e   1 × e x p r e s s i o n   l e v e l   o f   g e n e   1 + · · · · + β   g e n e   n × e x p r e s s i o n   l e v e l   o f   g e n e   n Here, ß is the regression coefficient in the multivariate Cox regression analysis. The patients in each cohort were divided into high- and low-risk groups, based on the median risk score of the training cohort. Univariate and multivariate Cox regression analyses were performed to evaluate the independent prognostic value of the risk signature. The R package “caret” was used to randomly split TCGA-PRAD into training and test cohorts at a 7:3 ratio. Statistical significance was set at p < 0.05.

The tumor mutation burden (TMB) of the TCGA-PRAD samples was visualized by the “maftools” R package (Mayakonda et al., 2018). TMB was defined as the number of base substitutions, deletions, insertions, and insertions across bases per megabase of the genome examined using non-synonymous and code-shifting indels under a 5% detection limit. Somatic alterations in PRAD driver genes were analyzed in samples with low- and high-risk scores.

The Gene_DE module of Tumor Immune Estimation Resource version 2 (TIMER2.0; http://timer.cistrome.org) (Li et al., 2020) was used to analyze the differentially expressed TIIC-related genes between the tumor and normal tissues in pan-cancer.

Seven methods, comprising XCELL, TIMER, QUANTISEQ, MCPcounter, EPIC, CIBERSORT, and CIBERSORT-ABS, were implemented to evaluate the extent of immune infiltration and its correlation to the risk score. The Estimation of Stromal and Immune cells in Malignant Tumors using the Expression Data (ESTIMATE) algorithm (Yoshihara et al., 2013) was used to predict tumor purity for each TCGA-PRAD sample.

The sensitivity of PRAD samples in the high- and low-risk-score groups to antitumor drug therapy was predicted by using the R package “pRRophetic” (Geeleher et al., 2014) to estimate the half-maximal inhibitory concentration (IC50) of each sample, based on the largest publicly attainable pharmacogenomics database: the Genomics of Drug Sensitivity in Cancer (GDSC) (www.cancerrxgene.org) (Yang et al., 2013) cell line expression spectrum.

To further explore the potential role of risk score in immunotherapeutic prediction, the Cancer Immunome Database (TCIA, https://tcia.at/home) provided comprehensive immunogenomic analyses of next-generation sequencing (NGS) data for 20 solid-tumor cancers from TCGA and other data sources. The immunophenoscore (IPS) was used as a novel and robust predictor of response to anti-cytotoxic T lymphocyte antigen-4 (anti-CTLA-4) and anti-programmed cell death protein 1 (anti-PD-1) antibodies (Charoentong et al., 2017) from the downloaded TCGA-PRAD datasets. The R package “ggpubr” was used to visualize IPS in the high- and low-risk groups. Furthermore, the expression levels of 47 immune checkpoint blockade-related genes in the high-and low-risk-groups were compared, and their correlations were visualized.

Molecular and functional relevance analyses of potential prognostic differentially expressed TIIC-related genes (PDEMRGs) were performed using Metascape (http://metascape.org) (Zhou et al., 2019). The activation of the hallmark pathway and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, described in the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) (Subramanian et al., 2005), was carried out to evaluate the enriched pathways in the high- and low-risk groups. To elucidate the functional annotation of each gene, including the risk signature, the gene set variation analysis (GSVA) R package was used to analyze the enrichment of the KEGG and gene ontology (GO) pathways.

R software (version 4.0.3) was used for all statistical analyses. The Wilcoxon test was used to compare the two groups. The R package “survivalROC” was used to calculate the area under the curve (AUC) using receiver operating characteristic (ROC) curves to identify the accuracy of the risk score. Kaplan–Meier curves with log-rank tests were used to compare survival rates. The chi-squared test was performed to correlate the risk-score subgroups with somatic mutation frequency, and Spearman analysis was used to compute the correlation coefficient. The “clusterProfiler”, “enrichplot”, “pheatmap”, and “ggplot2” R packages were used to visualize the results. Statistical significance was set at p < 0.05.

Five human prostate cancer cell lines (PC-3, DU-145, C4-2, 22RV-1, and LNCAP) were used to detect mRNA levels of activated dendritic cell-related genes. Among them, the PC-3 cells and DU-145 cells were purchased from American Type Culture Collection (ATCC), 22RV-1 cells and LNCAP cells were purchased from the Shanghai Fuheng Biotechnology Co., Ltd., and C4-2 cells were a gift from the Department of Endocrinology, the first hospital of Jilin University, Changchun, Jilin, China. All cell lines were cultured in Roswell Park Memorial Institute (RPMI-1640) medium supplemented with 10% fetal bovine serum (VivaCell, Shanghai XP Biomed Ltd.,). All media were supplemented with 5000 U/mL penicillin–streptomycin (Gibco). All cell lines were grown in a humidified atmosphere containing 5% CO₂ at 37°C. RNA was isolated using RNAiso Plus [Takara Biomedical Technology (Beijing) Co., China]. Total RNA (2.0 μg) was subjected to reverse transcription PCR (RT-PCR) using Hifair III Reverse Transcriptase (Yeasen, China) to obtain cDNA. cDNA was diluted 20-fold; then 6 μl was used for quantitative real-time polymerase chain reaction (qRT-PCR) using Hieff^® qPCR SYBR Green Master Mix (No Rox) (Yeasen, China). Gene expression levels were evaluated relative to GAPDH level and calculated using the 2^−ΔΔCt method. All samples were analyzed at least in triplicates. The primer sequences used for PCR were as follows: STX4, 5′- CGG​ACA​ATT​CGG​CAG​ACT​ATT -3′ (forward) and 5′- TTC​TGG​GGC​TCT​ATG​GCC​TT -3′ (reverse); UBE2S, 5′- CCG​ACA​CGT​ACT​GCT​GAC​C -3′ (forward) and 5′- GCC​GCA​TAC​TCC​TCG​TAG​TTC -3′ (reverse); TMEM93, 5′- GCC​GCC​GTC​CTG​GAT​TAT​T -3′ (forward) and 5′- GAG​GCG​AGC​AGG​TAG​AAG​AT -3′ (reverse); EMD, 5′- CCG​CCT​CCT​CTT​ATA​GCT​TCT -3′ (forward) and 5′- CTC​TGG​TAG​AGT​AAA​GCG​TCC​T -3′ (reverse); GCAT, 5′- CCT​CAG​CTC​TGT​CCG​CTT​TAT -3′ (forward) and 5′- GGA​TGC​CGT​CGA​TGA​TGG​AG -3′ (reverse); NUCB1, 5′- CAG​AAC​CAG​CAT​ACA​TTC​GAG​GC -3′ (forward) and 5′- AGT​GAC​TCC​AGA​TAA​CGC​CGT​C -3′ (reverse); and GAPDH, 5′- TCA​ACA​GCG​ACA​CCC​ACT​C-3′ (forward) and 5′- GCT​GTA​GCC​AAA​TTC​GTT​GTC-3′ (reverse).

We used the CIBERSORT algorithm to investigate 22 TIICs subsets in 550 samples (498 tumor and 52 normal) from the TCGA-PRAD dataset (Supplementary Table S1; Figure 1A). Ultimately, the abundance of 21 TIIC types from each sample was selected in each dataset (Supplementary Table S1; Figure 1A), excluding the naive T-cell type owing to its low abundance in all samples. A total of 78 samples under the threshold of the adjusted CIBERSORT p-value <0.05 were selected for subsequent analysis. The heatmap shows the TME patterns of 21 TIIC types in normal and tumor tissues (Figure 1B). Furthermore, the correlation matrix displays correlation coefficients between the 21 TIICs types, demonstrating a potential connection between these infiltrating immune cells in the TME (Figure 1C). Notably, CD8⁺ T-cells and regulatory T-cells (Tregs) had the strongest positive correlation (r = 0.51; p < 0.01), whereas follicular helper T-cells and resting memory CD4⁺ T-cells had the strongest negative correlation (r = −0.51; p < 0.01). In addition, the activated dendritic cells were positively correlated with several infiltrating immune cells, including memory B cells (r = 0.31; p < 0.01), naive B cells (r = 0.29; p < 0.01), and resting NK cells (r = 0.29; p < 0.01), suggesting their important roles in the TME (Supplementary Table S2).

We next established a WGCNA co-expression network through the “WGCNA” R package, selecting 19,560 genes from TCGA-PRAD patient samples after gene sequencing and preprocessing. A power of ß = 11 was used as the best soft-thresholding parameter, with an index of scale-free topologies of R ² = 0.90 (Figure 2A). Highly similar genes were assigned to the same module by the dynamic tree-cutting algorithm, and modules that met the hierarchical clustering analysis threshold of below 0.25 were clustered together. Consequently, 13 modules were identified in the resulting network and illustrated in a hierarchical clustering tree: MEmidnightblue, MElightyellow, MEred, MEblack, MEbrown, MEblue, MEgreenyellow, MEcyan, MEpurple, MEgrey60, MEmagenta, MEturquoise, and MEgrey (Supplementary Table S3; Figure 2B). Furthermore, the correlations between the 21 TIICs and each module were analyzed and presented as a heatmap (Figure 2C). Among these 13 modules, MEblue (r = −0.55, p = 1e⁻⁴⁰), MEgrey (r = −0.47, p = 1e⁻²⁸), and MEgreenyellow (r = −0.48, p = 2e⁻³⁰) modules were significantly negatively correlated with activated dendritic cells (Figure 2C). Eventually, we identified the MEblue module (a total of 1,143 genes) as having the strongest correlation coefficient among the analyzed results.

We selected the TCGA-PRAD dataset as the training cohort and the GSE116918 dataset as the validation cohort, the clinical features of which are listed in Table 1, and then extracted the expression levels of genes in the MEblue module from both TCGA-PRAD and GSE116918 datasets. Univariable Cox regression analysis of these genes revealed that 211 genes showed a significant prognostic value associated with PFS (p < 0.01) (Supplementary Table S4). LASSO Cox regression analysis of this group identified 12 genes with the highest coefficients (Figure 3A; Supplementary Table S5). Multivariate Cox regression analysis narrowed this group to six genes (STX4, UBE2S, EMC6, EMD, NUCB1, and GCAT) as the minimum set for constructing a TIIC-related genes risk signature (Table 2).

The risk formula was as follows: Risk score = (0.891 × expression level of STX4) + (0.618 × expression level of UBE2S) + (0.427 × expression level of EMC6) + (0.694 × expression level of EMD)—(0.628 × expression level of NUCB1)—(0.632 × expression level of GCAT). Among these genes, STX4 (syntaxin 4), UBE2S (ubiquitin conjugating enzyme E2 S), EMC6 (ER membrane protein complex subunit 6 or transmembrane protein 93, TMEM93), and EMD (emerin) were identified as high-risk genes, whereas NUCB1 (nucleobindin 1) and GCAT (glycine C-acetyltransferase) as low-risk genes. Survival analysis showed that TCGA-PRAD patients with higher expression levels of high-risk genes had poorer prognoses than those with lower expression levels (Supplementary Figures S1A–D). The PFS of patients with higher expression levels of low-risk genes was greater than those with lower expression levels (Supplementary Figures S1E–F).

Furthermore, we analyzed the mRNA expression levels of the six TIIC-related genes from TCGA pan-cancers in the TIMER database, and found that they were significantly increased in tumor tissues compared to the adjacent normal tissues found in various cancers (p < 0.05) (Supplementary Figures S2, S3). For example, the high-risk genes EMC6 and UBE2S were expressed to a greater extent in TCGA-PRAD tumor tissues than in normal tissues. Although STX4 was not overexpressed in TCGA-PRAD, its mRNA level was higher in stage IV than in stages I–III (Supplementary Figure S1G). Intriguingly, we found that the mRNA levels of low-risk genes NUCB1 and GCAT were elevated in tumor tissues more than those in normal tissues, but they were lower in stage IV than in stages I–III (Supplementary Figures S1H–I). Together, these results suggest that the three genes may play essential roles in the advanced disease stage and thus validate the combined influence of the six activated dendritic cell-related genes on prostate cancer risk. Presumably, EMC6 and UBE2S are overexpressed in tumor tissues compared with normal tissues. They are thus likely to be oncogenic from the start, whereas STX4 probably do not play a role during the initial stages of cancer, and NUCB1 and GCAT are the significant players during this phase.

To further assess outcome prediction, we calculated the risk scores for each patient in the training cohort using the six-gene formula model and divided them into high- and low-risk groups based on the median risk score cutoff value of 0.9773 (p < 0.001). Kaplan–Meier analysis showed that PFS was shorter in the high-risk group than in the low-risk group in both the training (p < 0.001) and validation (p = 0.031) cohorts (Figures 3B, C). Scatter plots and risk curves indicated that the signature gene risk score was negatively proportional to the PFS of PCa patients (Figures 3D, E). These results suggest that the risk signature of the six TIIC-related genes has an excellent prediction performance and thus demonstrates the potential as a prognostic factor in PCa patients.

Next, a 5-year ROC curve analysis was performed, revealing that the ROC of the risk score (5-year PFS AUC, 0.753) was significantly higher than that of the prognostic-related clinical parameters, such as stage (0.569) (Figure 3F), validating the accuracy of the risk model. Furthermore, heatmap analysis showed that, in the high-risk group, the mRNA expression levels of STX4, UBE2S, EMC6, and EMD were higher, whereas those of NUCB1 and GCAT were lower than those in the low-risk group (Figures 3G, H), suggesting their prognostic-specific roles in PCa patients. In addition, univariate and multivariate regression analyses indicated that the risk signature was an independent prognostic predictor for PFS, with hazard ratios of 1.076 (95% CI:1.053–1.100) in the univariate analysis and 1.068 (95% CI:1.044–1.093) in the multivariate analysis (Figures 3I, J). To further validate the stability of the model, we randomly split TCGA-PRAD into training and test cohorts at a 7:3 ratio and constructed another prognostic model with different variables. We found that the new model had four genes similar to our previous risk signature and accurately predicted the prognosis in the TCGA test and GSE116918 cohorts (Supplementary Table S6; Supplementary Figures S4A–F).

In addition, we examined the mRNA expression levels of the six signature genes in five different PCa cell lines and found that UBE2S had the highest expression level among the six genes (Figure 4).

Next, we used a heatmap to visualize the distribution difference of clinical variables between the high- and low-risk groups and found that most patients in the high-risk group were older (>65 years), T-Stage scores (≥3), and N1-stage and Gleason scores (≥7) than those in the low-risk group (Figure 5A). Bar plots also confirmed the proportion difference in clinical subtypes based on age and T-stage, N-stage, and Gleason scores in the high- and low-risk groups (Figures 5B–E). Furthermore, we performed subgroup analysis to determine whether our risk signature could identify different prognoses. When patients were classified based on age, our risk signature accurately predicted patient outcomes, with higher scores indicating poorer outcomes (Supplementary Figures S5A, B). The risk signature was consistently capable of prognostically predicting patients in the T3 category (Supplementary Figure S5D), those with N0 status (Supplementary Figure S5F), and those with a Gleason score of 7 (Supplementary Figure S5I). Notably, the insignificant prognostic prediction of the risk signature in some of the clinicopathological parameter subgroups (i.e., Gleason score of 6, T-stage of T4, and N-stage of N1) was likely due to the relatively low number of cases in the study. However, it is worth noting that our risk signature demonstrated better predictions for PRAD patients in the late stages of PCa than for those in the early stages (i.e., Gleason score >6 and > T2 stage). These findings, combined with the results of the univariate and multivariate regression analyses (Figures 3I, J), indicate the statistical and clinical significance of the risk signature as a prognostic, predictive indicator.

PCa frequently exhibits genomic alterations (to the AR axis, ETS family, TP53, PTEN, or RB) (Cancer Genome Atlas Research, 2015; Abida et al., 2019). These include alterations in genes involved in biochemical pathways (the PI3K/AKT/MAPK pathway and cell cycle-related pathways) (Cancer Genome Atlas Research, 2015; Abida et al., 2019; Powles et al., 2022), epigenetic changes (Cancer Genome Atlas Research, 2015; Abida et al., 2019), alterations in DNA repair pathways, including homologous recombination repair (HRR) and mismatch repair (MMR) (Cancer Genome Atlas Research, 2015; Abida et al., 2019; Devlies et al., 2020), and single-nucleotide variants (SNVs) (AR, TP53, PI3KCA, BRCA2, PTEN, APC, CDK12, and ATM) (Robinson et al., 2015; Abida et al., 2019; Powles et al., 2022). TMB is calculated based on the somatic mutation frequency and has been proposed as a predictor of immunotherapy efficacy in various cancers, including bladder cancer, NSCLC, small cell lung cancer, and melanoma (Rizvi et al., 2015; Yarchoan et al., 2017). We compared the TMB of patients in the low- and high-risk groups and found that the TMB of the high-risk group was significantly higher than that of the low-risk group (p = 1.4e⁻⁶) (Figure 6A). The risk score positively correlated with TMB (R = 0.25, p = 6.5e⁻⁸; Figure 6B).

Furthermore, we examined somatic variations in driver genes in PCa between the low-risk and high-risk groups, and the top 20 genes with the highest mutation frequencies were analyzed (Figures 6C, D). We found that the SPOP (speckle-type POZ protein) gene was the most significantly mutated gene (SMG) in both risk-score subgroups (11% and 9%) (Figures 6C, D), confirming previous reports that SPOP is the most frequently mutated gene in PCa (Jin et al., 2021). TP53 (12% vs. 7%) had higher somatic mutation rates in the high-risk group, whereas RYR2 (4% vs. 1%) had higher somatic mutation rates in the low-risk group. These findings demonstrate the genetic profiles underlying the intrinsic connection between activated dendritic cell infiltration and somatic variants in PCa immunotherapy.

Then, the patients were separated into a high-TMB group (n = 231) and a low-TMB group (n = 241) based on the optimal cutoff value of TMB (cutoff value = 0.474). Kaplan–Meier curves showed that patients in the low-TMB group had better PFS than those in the high-TMB group (log-rank test, p ≤ 0.001; Figure 6E). To evaluate the synergistic effect of the TMB grouping and risk-score grouping in the prognostic stratification, patients were divided into subgroups of high-TMB and high-risk, high-TMB and low-risk, low-TMB and high-risk, and low-TMB and low-risk groups based on the optimal cutoff value of TMB and median risk score cutoff value. The TMB status did not affect the survival prognosis prediction based on the risk-score group. However, the risk-score subgroup showed significant survival differences in both the low- and high-TMB groups (log-rank test, high-TMB and high-risk vs. high-TMB and low-risk, p < 0.001; low-TMB and high-risk vs. low-TMB and low-risk, p < 0.001; Figure 6F). Notably, the low-TMB and low-risk subgroups had the best PFS rates, whereas the low-TMB and high-risk subgroups had the worst PFS rates (Figure 6F).

To investigate the intrinsic and intimate connection between the risk signature and TIICs in the TME, we performed a correlation analysis of the risk score in the TME context of PRAD and found that the risk score had a strong, negative correlation with subpopulations of naive B cells (CIBERSORT, R = −0.2, p < 0.01), plasma B cells (CIBERSORT, R = −0.4, p < 0.01), and neutrophil cells (QUANTISEQ, R = −0.35, p < 0.01) but was positively correlated with an abundance of B cells (XCELL, R = 0.12, p < 0.01), memory B cells (CIBERSORT-ABS, R = 0.22, p < 0.01), M1 Macrophages (XCELL, R = 0.2, p < 0.01), activated NK cells (CIBERSORT-ABS, R = 0.17, p < 0.01), activated myeloid dendritic cells (CIBERSORT-ABS, R = 0.19, p < 0.01), activated myeloid dendritic cells (XCELL, R = 0.2, p < 0.01), and follicular helper T-cells (CIBERSORT, R = 0.19, p < 0.01) (Supplementary Table S7; Supplementary Figures S5A, B). We also analyzed the correlation between the risk signature and immune infiltration (Figure 7A). ESTIMATE analysis showed that the immune (p < 0.01), stromal score (p < 0.05), and ESTIMATE scores (p < 0.01) were significantly higher in the high-risk group than in the low-risk group (p < 0.01) (Figure 7B).

In addition, we determined the optimal cut-off value of the ESTIMATE score (−390.9762) using the minimum p-value method and classified the patients into a high-ESTIMATE group (n = 175) and low-ESTIMATE group (n = 318). Kaplan–Meier curves showed that patients in the low-ESTIMATE group had better PFS rates than those in the high-ESTIMATE score group (p = 0.011; Figure 7C). Furthermore, we analyzed the synergistic effect of the ESTIMATE score and risk-score grouping in the prognostic stratification. Kaplan–Meier curves indicated that the ESTIMATE score did not affect the survival prognosis prediction based on the risk-score subgroup. The risk-score subgroup demonstrated significant survival differences in low and high ESTIMATE subgroups (log-rank test, high-ESTIMATE and high-risk vs. high-ESTIMATE and low-risk, p < 0.001; low-ESTIMATE and high-risk vs. low-ESTIMATE and low-risk, p < 0.001; Figure 7D). Notably, the low-ESTIMATE and low-risk subgroups had the best PFS rates, and the high-ESTIMATE and high-risk subgroups had the worst PFS rates (Figure 7D).

We further analyzed the association between the risk score and efficacy of antitumor therapeutic drugs for PCa. Although many immune checkpoint blockade-related genes (e.g., PDCD1 and CTLA4) showed a significantly positive correlation with the risk score (Figure 8A), the scores of the IPS-PD1 blocker, IPS-CTLA4 blocker did not reveal significant differences between the high-risk and low-risk groups (Supplementary Figures S7A–C). Intriguingly, the IPS score (PD-1 negative and CTLA-4 negative; Figure 8B) of the high-risk group was higher than that of the low-risk group, suggesting that patients in the high-risk group could have benefited from immune checkpoint blockade (ICB) treatment instead of PD1/CTLA4 immunotherapy. Furthermore, we used the pRRophetic package to evaluate the sensitivity of the risk score to antitumor therapeutic drugs. We found that three common PCa drugs (bicalutamide, docetaxel, and rucaparib/AG.014699) showed different sensitivities in the high-risk and low-risk groups. Furthermore, patients in the high-risk group had higher half-inhibitory concentration (IC50) values than those in the low-risk group, indicating that they were less sensitive to these three antitumor therapeutic drugs (Figures 8C–E).

We performed gene set variation analysis (GSVA) in the training cohort to analyze the signal pathways activated in high-risk or low-risk groups (Figures 9A, B). In the high-risk group, E2F-regulated, DNA repair, MYC-regulated, UV-activated, glycolysis, and p53-mediated signaling pathways were activated. These signaling pathways are involved in advanced disease progression (Mandigo et al., 2022), sensitivity to drug therapies (Wei et al., 2021), drug resistance (Jividen et al., 2018), tumor immune infiltration (Elliott et al., 2019), tumor proliferation (Elliott et al., 2019), and immune resistance (Cascone et al., 2018). In the low-risk subgroup, the WNT, PPAR, protein secretion, cholesterol homeostasis, and androgen response signaling pathways were elevated. These signaling pathways are involved in the regulation of the disease metastasis (Leibold et al., 2020), the neuroendocrine differentiation (Liu et al., 2019), and the cell cycle, proliferation, and migration (Aurilio et al., 2020). These results show differences in the biological processes between the high-risk and low-risk groups.

The endomembrane system organization process was the most significantly enriched gene ontology biological process (GO-BP) among the six TIIC-related genes (Supplementary Figure S5A). We also used gene set enrichment analysis (GSEA) to analyze the functional enrichment activated by the six TIIC-related genes and identified significant enrichments in both the GO pathways (Supplementary Table S8; Supplementary Figures S5B–G) and KEGG pathways (Supplementary Table S9; Supplementary Figures S6A–F).