We downloaded the transcriptome data (read counts and transcripts per million (TPM)) and corresponding clinical data of PCa from the TCGA database using the R package “TCGAbiolinks”, and randomly divided the tumor data (n = 500) into five parts using R package “Caret” and setting seed as “1”, one of which was the internal validation set (n = 100) and the rest (n = 400) the training set. The external validation set EGAS00001002923 (Gerhauser et al., 2018) with RNA-seq expression (reads per kilobase per million mapped reads (RPKM)) and corresponding clinical data were obtained through the cBioPortal (Prostate Cancer (DKFZ, Cancer Cell 2018)). Data had to meet the following inclusion criteria at the same time: the availability of RNA-seq sequencing data, clinical progression-free survival data and BCR-related indicators. In the aggregate, 552 TCGA-PRAD samples across 52 normal tissues and 500 prostate cancer tissues, of which 400 were a training set and 100 were an internal validation set, as well as 82 samples from EGAS00001002923 used as an external validation set, were collected for analysis. Furthermore, the gene signatures of 28 tumor-infiltrating lymphocytes were captured from the TISIDB database (http://cis.hku.hk/TISIDB/), listed in Supplementary Table S1. A total of 797 LRGs from 23 lysosome-related pathways were obtained from the Molecular Signatures Database, listed in Supplementary Table S2.

To compare the composition and differences among 28 kinds of tumor-infiltrating immune cells between cancer samples and normal samples from TCGA-PRAD, we quantified the ssGSEA scores by using an ssGSEA algorithm and R package “GSVA”. Of note, the ssGSEA scores were normalized to a 0–1 interval and processed by the mean value. RNA expression was expressed as TPM.

PCa samples from the TCGA database were divided into high-immune and low-immune groups based on the median ssGSEA score. Next, R package “DESeq2” was used to screen out differentially expressed genes (DEGs) between two immune groups; 1,093 IRGs with the adjusted p-value < 0.05 and the absolute value of Log2FoldChange >1 was considered as statistically significant and were visualized by a volcano plot. After that, a gene ontology (GO) enrichment analysis and its visualization were applied for IRGs by using R packages “topGO” and “clusterProfiler”. RNA expression was expressed as read counts.

We generated a Venn diagram to show LIRGs, which were categorized as the shared gene between IRGs and LRGs. Next, we performed a heat map to compare levels of LIRGs between low and high-immune groups. In order to explore the correlation of individual LIRGs and PCa prognosis, the patient progression-free interval (PFI) was determined by the Kaplan–Meier (KM) method and Log-rank test between the high- and low-expression groups based on the scanning cut points with R package “survival” and “survminer”.

A univariate Cox regression analysis was used to screen out progression-associated factors among the training set and was shown in a forest plot using R package “forestplot”. Moreover, we used the least absolute shrinkage and selector operation (LASSO) algorithm to further narrow down candidate factors. Then, a multivariate Cox regression analysis was performed to select the factor with independent prognostic value and determine the coefficients of a linear risk score formula. The linear risk score formula of the whole PCa can be calculated as follows: R i s k   s c o r e = ∑ i = 1 n X i × β i

In this equation, n represents the number of prognostic factors, X represents the expression of genes or the Gleason score and ß represents the regression coefficients of prognostic factors. The proportional hazards (PH) assumption can be checked using R package “survival”. Meanwhile, the training set was grouped into high and low-risk groups based on the cut point of the risk score using R package “timeROC”. A KM curve by log-rank test and a time-dependent receiver operating characteristic (ROC) curve were performed to evaluate the prognostic ability and performance of the prognostic model. According to the results of the multivariate analysis, we applied R package “rms” to establish a nomogram for visualizing the model, which contributes to the guideline of clinical decision-making. The area under the curve (AUC) by R package “timeROC”, the concordance index (C-index) and the calibration curve by R package “rms” were used to assess the model’s predictive accuracy. Furthermore, the performance of the model was verified in the internal and external validation set using the same methods. In the external validation set, we used biochemical recurrence (BCR) to represent progression. The RNA-seq expression mentioned above was processed by Log₂ (TPM+1).

All statistical analyses were performed under R environment (R version 4.2.1). The unpaired two-sample Wilcoxon test (also known as the Mann–Whitney test) was utilized to compare two independent groups. Comparison of categorical variables was evaluated by Fisher’s exact test. R packages “ggplot2”, “ggpubr” and “pheatmap” were applied for different plots. p < 0.05 was set as statistically significant in our study.

We obtained data from 500 prostate cancer samples and 52 normal samples from the TCGA-PRAD cohort. Demographics, clinical characteristics and progression status are presented in Table 1. Compared to the non-progression group, mean age is higher in progression group. Besides, the Gleason score shows the most significant difference between the two groups by univariate comparisons of the demographic and clinical characteristics. To assess the unique immune cell infiltration of prostate cancer, which is distinct from normal tissues, we applied the ssGSEA algorithm based on the RNA-sequencing data. The ssGSEA score of 28 immune-related cells in both normal and tumor samples is presented in Figure 1A. The results show that the infiltration of most immune cells (20/28, 71.4%), especially most subtypes of T cells, macrophages, mast cells, plasmacytoid dendritic cells, monocytes and myeloid-derived suppressor cells (MDSC), was significantly different between the normal and tumor samples (p < 0.05), which had high potential relativity between the immune microenvironment and prostate cancer. Additionally, the boxplot demonstrated that prostate cancer samples had a lower ssGSEA score, indicating a suppressed immune environment in prostate cancer (Figure 1B). Prostate cancer samples (n = 500) were thereafter divided into two immune subgroups based on the median ssGSEA score; the DEGs are listed in Supplementary Table S3. Compared to the high-immune group, the mRNA expression of most genes in the low-immune group was downregulated (Figure 1C). Of note, these differentially expressed genes in two immune subgroups were considered as IRGs (n = 1,093). The GO enrichment analysis of 1,093 IRGs is shown in Figure 1D. After clustering the GO cellular component enriched pathways, some could be classified into the primary lysosome azurophil cytoplasmic, indicating that some of the IRG products are located in lysosomes or are related to lysosomes when they perform their function.

As shown in the Venn diagram, a total of 49 (49/1841, 3%) shared genes between IRGs and LRGs were screened out and regarded as LIRGs (Figure 2A). Figure 2B demonstrates the 49 LIRG expression patterns between two immune subgroups. To assess the prognostic relevance of the 49 candidate LIRGs, we utilized the KM method and found that more than half of the LIRGs (29/49, 59.2%) were associated with the prognosis (all p < 0.05) as shown in Figure 2C. As we can see, the high expression of most LIRGs (23/29) predicted a poor prognosis.

A univariate Cox regression analysis was performed with the Gleason score and the expression profiles of 49 candidate LIRGs in the training set (Figure 3A). In addition to the Gleason score, 15 candidate LIRGs were found to have a significant correlation with progression free interval (all p < 0.05, marked in red); the baseline is shown in Supplementary Table S4. Of note, these 16 prognostic signatures were risk factors (HR > 1). To further narrow down candidate prognosis signatures, three candidates, including Gleason score (coefficient = 0.6509428), NCF1 (coefficient = 0.3190647) and IFI30 (coefficient = 0.4826147), were identified based on the optimum λ value 0.02010609 using a LASSO regression analysis (Figures 3B, C). Then, a multivariate Cox regression analysis further screened out the prognosis-associated signature and established a linear risk score formula. As shown in Figure 3D, the risk score was calculated as follows: R i s k   s c o r e = 0.7470 × G l e a s o n   s c o r e + 0.4020 × Log 2 T P M N C F 1 + 1 + 0.6365 × Log 2 ( T P M I F I 30 + 1

To evaluate the performance of the predictive model constructed by using the risk score, four time-dependent ROC curves were drawn, where the AUC values at 1, 3, 5 and 10 years were 0.790, 0.779, 0.779 and 0.894, respectively (Figure 3E), indicating excellent predictive ability. Patients with prostate cancer were subsequently divided into high and low-risk groups based on the risk score of 6.275728; the cutoff value of the risk score was determined by the ROC analysis (Figure 3F).

As shown in the box plot, Gleason scores increased and the expression of NCF1 and IFI30 were upregulated, along with the increasing risk score (Figure 4A). On the basis of the progression free interval analysis, patients in the high-risk group were substantiated to have poorer prognoses by the log-rank test (p < 0.0001; Figure 4B). Likewise, in the cumulative hazard analysis, patients with higher risks showed a higher cumulative hazard (Figure 4C). Meanwhile, in order to demonstrate the evaluation of patients’ PFI probability, the nomogram model was established by combining the three prognosis-associated features mentioned above (Figure 4D). Moreover, both 3- and 5-year PFI predictions by nomogram were highly consistent with the actual observation of patients with prostate cancer, as confirmed in the calibration plots (Figure 4E). These results verify that this model, constructed by three prognosis-associated features, is predictive of poor prognoses. Interestingly, the AUC of this risk model was greater than that of the model for which only the Gleason score was considered, whether at 1, 3, 5 or 10 years (Figures 4F, G), suggesting that NCF1 and IFI30 can refine the abilities of prognostic prediction by using the Gleason score alone, which has been the most well-accepted and strongest prognostic predictive tool in prostate cancer.

As mentioned in the Methods section, except for the training set, the remaining 100 cases from TCGA-PRAD were used as the internal validation set (Supplementary Table S5), while 82 cases from EGAS00001002923 were used as the external validation set (Supplementary Table S6). In addition, we also considered the entire TCGA-PRAD database (n = 500) as a kind of data cohort for the internal validation set (Supplementary Table S7). First, the calibration curves based on three validation sets showed good agreement between the actual and predicted 3- and 5-year PFI, which illustrated that the risk model was an excellent predictor (Figures 5A–C). Further analysis of AUC by ROC curves ranged from 0.658 to 0.812 (Figures 5D–F), and proved the robust predictive capacity of our model. Meanwhile, a univariate regression analysis was also used to verify the prognostic role of risk score in prostate cancer (Figures 5G–I), revealing that patients with a higher risk score had a poor prognosis.