Establishment of a novel prognostic prediction model through bioinformatics analysis for prostate cancer based on ferroptosis-related genes and its application in immune cell infiltration.

Background: Ferroptosis-related genes (FRGs) play vital roles in survival and prognosis of prostate cancer (PCa) patients. We establish a ferroptosis-related prediction model through bioinformatics analysis for overall survival (OS) and disease-free survival (DFS), so as to evaluate the clinical survival status through the characteristics of immune cell infiltration (ICI), which could provide information for treatment monitoring.
Methods: At first, 268 FRGs were obtained from previous studies. Differentially expressed FRGs were identified based on The Cancer Genome Atlas (TCGA) database, and FRG enrichment analysis was performed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). We then performed univariate, least absolute shrinkage and selection operator (LASSO), and multivariate Cox regression analyses to establish OS- and DFS-related prognostic prediction models. The association of the model and clinicopathological features was further analyzed. Subsequently, unique genomic signatures of immune cell subsets were obtained through the KEGG database. Based on specific genes associated with ferroptosis and their association with ICI, immune infiltration was assessed in patients in different risk groups.
Results: We constructed an OS- and an DFS-prognostic model through bioinformatics analysis. The predicted values of OS and DFS-related models were higher in T3-4 than in T1-2 (P=0.0057, P<0.001), and the predicted value of the DFS model in N0 stage was higher than that in N1 stage (P=0.0136). Results of Single-sample gene set enrichment analysis (ssGSEA) on the basis of the KEGG dataset showed p53 signaling being the most enriched signal in the high-risk group, while endocytosis was the most enriched signal in the low-risk group. M2 macrophages (P=0.007) and neutrophils (P=0.024) were enriched in the high-risk group, and CD4-activated memory T cells were significantly accumulated in the low-risk group (P=0.017). Conclusions: The OS- and DFS-related model based on FRGs and ICI create new insights into the disease state assessment of PCa patients., which may aid in the development of individualized and precise treatment in the future. 2022 Translational Andrology and Urology. All rights reserved.

Introduction

Prostate cancer (PCa) is one of the most prevalent male malignancies worldwide according to the latest cancer statistics for 2022 (1). For patients with localized PCa, 20–25% of patients may still experience biochemical recurrence (BCR) during follow-up despite surgery or other treatments, and severe cases may eventually progress to castration-resistant prostate cancer (CRPC) with increasing specific mortality and survival risk (2). Therefore, the development of new prognostic identification markers is crucial for monitoring the clinical survival status of PCa patients, especially high-risk PCa patients.

In 2012, in order to describe the morphology of cell death induced by erastin, a small molecule binds and inhibits voltage-dependent anion channels, Dixon et al. coined the term ferroptosis. Ferroptosis is iron ion-catalyzed necrotic cell death induced by inhibiting cysteine import, resulting in glutathione depletion and inhibition of glutathione peroxidase 4 (GPX4) (3). In recent years, studies have reported that ferroptosis regulates the progression of various cancers by inhibiting tumor cell proliferation. Whether the process is genetically determined or pathologically triggered remains unclear. Under the catalysis of divalent iron or ester oxygenase, the highly expressed unsaturated fatty acids on the cell membrane undergo lipid peroxidation and lipid-reactive oxygen accumulation, leading to cell ferroptosis (4). Currently, the role of ferroptosis in tumor development, treatment, and prognosis has become an active area of research. PCa is considered one of the diseases associated with ferroptosis (5). Dyslipidemia was demonstrated to be the most commonly acquired metabolic phenotype in androgen receptor inhibitor (ARI)-resistant cells. In the treatment of CRPC, adjusting the ratio of saturated to unsaturated phospholipids increases the sensitivity of CRPC cells to ferroptosis and increases cytotoxicity (6). Erastin and the GPX4 inhibitor RSL3, 2 ferroptosis inducers (FINs), have been shown to slow tumor cell growth in vivo and significantly reduces tumor cell migration in vitro, suggesting that treatment-resistant PCa cells are highly sensitive to these 2 iron death inducers (5). Bordini et al. confirmed that androgen receptor-positive PCa cell lines VCaP and LNCaP have high iron sensitivity and eventually cause cell membrane protein damage, suggesting that ferroptosis is associated with poor prognosis in prostate tumors. Furthermore, in vitro experiments confirmed that bicalutamide-iron conjugate resulted in protein oxidation and successfully inhibited tumor expansion, whereas the single compound was ineffective (7). The above studies have confirmed that ferroptosis may become an entry point for the prognosis of PCa survival.

In the past decade, cancer immunotherapy research has made historic breakthroughs (8). Tumor microorganisms act on the tumor-associated immune microenvironment (TAIM) by inducing immunogenic cell death, assisting pattern recognition receptor-mediated signaling pathways, and bacterial antigens mimicking tumor antigens, which are of great significance for tumor development and immunotherapy response (9). A suppressive TAIM, characterized by immune cell and stromal cell infiltration, severely affects tumor proliferation, metastasis, recurrence, and immunotherapy resistance, as seen in studies of a variety of cancers (10). Based on the TCGA database, a variety of different immune cells can be detected from PCa tissue (11). So far, the pattern of immune cell infiltration (ICI) and clinical outcome in PCa are poorly characterized. An in-depth understanding of the TAIM is urgently needed to identify the PCa subtypes that contribute to differential prognosis and lay the groundwork for improved immunotherapy.

The above evidence shows that the expression of ferroptosis-related genes (FRGs) and ICI are closely related to the local growth and metastasis of PCa, which may has certain significance for predicting the individualized health prognosis of patients. For the first time, this study evaluated the individual status of PCa patients from multiple perspectives by constructing a prognostic characteristic model of overall survival (OS) and disease-free survival (DFS) with the combination of analysing the immune infiltration landscape analysis, in order to supplement the existing clinical staging system and promote the precision of cancer treatment. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-22-454/rc).

Methods

Based on bioinformatics analysis, this study firstly searched for FRG templates and screened out the genes in the group. Through univariate, LASSO, and multivariate regression analysis, a predictive model was constructed, and then the model was varified by means of clinicopathological samples; then the application of FRGs in immune infiltration was tested.

Data acquisition

A total of 268 FRGs were extracted from previous study (12). After intersecting these genes with the expression matrix of PCa transcriptome data in The Cancer Genome Atlas (TCGA; https://tcga-data.nci.nih.gov/tcga/), a total of 249 genes were involved in the follow-up analysis.

PCa patient cohort data (ReadCounts) were downloaded from TCGA database (https://portal.gdc.cancer.gov/) and be served as a training set. The RNA-Seq expression matrix of 500 patients was obtained using the R package TCGAbiolinks. For the purpose of obtaining comparable standardized databases, fragments per kilobase of transcript per million fragments mapped (FPKM) values of RNA-Seq were log2 transformed. Finally, 498 PCa tissues and 52 adjacent normal tissues were obtained. The clinical information of 500 patients with PCa was collected (). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Table 1

Major demographic and clinicopathological characteristics of PRAD cases (n=500)

Patients characteristic	Number (%)
Age (years)
≥65	167 (33.4)
<65	333 (66.6)
T stage
T1/T2	188 (38.1)
T3/T4	305 (61.9)
N stage
N0	348 (81.5)
N1	79 (18.5)
Gleason score
≥7	455 (91.0)
<7	45 (9.0)

PRAD, prostate adenocarcinoma.

Discernment of differentially expressed fRGs between PCa and adjacent non-tumor tissue

Differential expression analysis of FRGs in PCa tissues and adjacent non-tumor tissues was performed utilizing the Limma package in R, with thresholds of |log2 fold change (FC)| >1. FDR <0.05 was considered statistically significant. The results were presented in 3 forms: heat map, volcano plot, and boxplot.

Enrichment analysis

Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was forthputte with R package clusterProfiler. The results were represented by dot plots.

Construction and validation of a prediction model based on FRGs

To identify OS- and DFS-related FRGs in PCa, we firstly executed univariate Cox analysis to initially screen for statistically significant genes. Least absolute shrinkage and selection operator (LASSO) was processed so as not to overfit and remove tightly correlated genes. Subsequently, we used multivariate Cox regression with the R package Glmnet to build the OS- and DFS-related prediction model. The formula of the risk score was as follows:

where Coefm represents the risk coefficient and Expm represents the relative mRNA expression of each FRG. We stratified patients with PCa into high-risk and low-risk group based on median risk scores obtained from database data. Survival status was assessed by the Kaplan-Meier (K-M) method. Sensitivity and specificity of survival predictions were tested by receiver operating characteristic (ROC) analysis using the R package survivalROC. The area under the ROC curve (AUC) was calculated as an indicator of the accuracy of the predictions for both groups. The K-M curves of each single gene in the OS and DFS prediction model were analyzed separately to assess their correlation. Finally, we compared the FRGs in the OS and DFS prediction model with the clinicopathological parameters of PCa patients.

Relationships of FRGs with ICI

Single-sample gene set enrichment analysis (ssGSEA) was used to quantify the relative infiltration of immune cell types in the tumor microenvironment. The unique genomes of each independent immune cell subsets were obtained from the KEGG (https://www.genome.jp/kegg/) database. Gene set variation analysis (GSVA) was implemented with the help of R package GSVA to assess the relative abundance of immune cells between different risk groups. To standardize the data, 1 was set as the maximum score and 0 was set as the minimum score.

The CIBERSORT (https://cibersort.stanford.edu/) algorithm was used to assess the richness of ICI in different groups of patients. The Estimate algorithm was used to obtain the difference of Stromasignature and ImmuneSignature between high-risk group and low-risk group.

Statistical analysis

All data analysis was performed using R4.1.1 (https://cran.r-project.org). Survival curves for OS- and DFS-prognostic analysis were produced by the K-M method, and the log-rank test was used to acknowledge the distinguishment of differences. Univariate and multivariate Cox regression models were used to calculate hazard ratios (HR), and forest plots visualized the coefficients of these regression models. LASSO regression was applyed to avoid overfitting during the construction of the final prediction model. Continuous variables were compared between the 2 groups by the Wilcoxon rank test. Comparisons of categorical variables between 2 groups were performed by the chi-square test. The Kruskal-Wallis test was used to compare differences among 3 or more groups. As for ROC, models with AUC greater than 0.75 perform better. Statistical significance was defined as |log2 FC >1.0| and adjusted FDR P value <0.05.

Results

Differentially expressed FRGs between PCa and adjacent non-tumor tissues

A total of 500 patients with primary localized PCa containing complete RNA-Seq transcript data and clinical follow-up information including age, tumor stage, and gleason score were enrolled in this study. Among the 249 FRGs included in the study, 31 were differentially expressed ().

Figure 1

Differential expression of ferroptosis-related genes in prostate cancer and adjacent tissues. (A) Heat map of the differential expression of FRGs in 500 patients, where red represents positive correlation and blue represents negative correlation. (B) Volcano plot of 249 FRGs of TCGA-PRAD. Red indicates high expression, blue indicates low expression, and gray indicates that the difference is not statistically significant. (C) Expression patterns of differentially expressed FRGs in PCa and paired non-tumor samples. Red and blue indicate tumor tissue and adjacent normal tissue, respectively. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. FDR, false discovery rate; FRGs, ferroptosis-related genes; TCGA-PRAD, The Cancer Genome Atlas-Prostate Adenocarcinoma; PCa, prostate cancer.

Volcano plots and boxplots were generated to demonstrate the expression patterns of FRGs differentially expressed between PCa and adjacent normal tissues ().

GO and KEGG enrichment analysis of differentially expressed FRGs

GO (http://geneontology.org) and KEGG (https://www.genome.jp/kegg/) enrichment analysis was performed based on differentially expressed FRGs. According to the results of GeneRatio, most of the GO terms were enrich for several metabolism processes, such as reactive oxygen species metabolic process, fatty acid metabolic process, and response to oxidative stress. Interestingly, the KEGG terms were closely associated with microRNAs in cancer, ferroptosis, arachidonic acid metabolism, and thyroid synthesis hormone. Dot plots were used to show the relationship between FRGs and GO and KEGG enrichment analysis ().

Figure 2

GO enrichment analysis (A) and KEGG enrichment analysis (B) of the differential expression of FRGs. BP, biological processes; CC, cellular components; MF, molecular functions; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FRGs, ferroptosis-related genes.

Establishment and validation of OS and DFS models in PCa patients based on FRGs

In the initial anlaysis by univariate Cox regression, a total of 31 FRGs were significantly correlated with OS, and 6 of them (AKR1C1, MT3, SLC2A6, AURKA, AKR1C3, TFRC) were shown in the forest plot (). Hub genes were selected using LASSO to screen the colinearity and reduce the risk of overfitting simultaneously. The LASSO Cox regression method constructed a 21-mRNA signature ().

Figure 3

OS-related prognostic model of prostate cancer patients. (A) Univariate Cox regression analysis of 31 ferroptosis-related genes associated with OS in PCa patients. (B,C) A 22 mRNA signature was constructed using the LASSO Cox regression method. (D) Kaplan-Meier plot shows that the OS of patients in the high-risk group is significantly lower than that of the patients in the low-risk group. (E) ROC subject curve: the curve of AUC at 1 year, 3 years, and 5 years, where AUC >0.75 predicts the results with statistical significance. (F,G) OS-related prediction model distribution of patients in TCGA dataset. (H) Heat map showing the correlation of 4 genes between the high- and low-risk groups. AUC, area under the ROC curve; OS, overall survival; PCa, prostate cancer; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

In addition, 4 genes, ALOX12, AKR1C2, SRC, and BNIP3, were identified as predictors and confirmed in the multivariate Cox regression analysis to establish a prediction model: OS-related prediction model = (2.5213096 × expression value of ALOX12 + (0.9157942 × expression value of AKR1C2) + (1.106994 × expression value of SRC) + (0.7573778 × expression value of BNIP3). Among them, 2.5213096, 0.9157942, 1.106994 and 0.7573778 are the distribution coefficients of the predictors ALOX12, AKR1C2, SRC and BNIP3 in the model, respectively.

Subsequently, according to their median risk score, we divided all the cases in The Cancer Genome Atlas-Prostate Adenocarcinoma (TCGA-PRAD) cohort into high- and low-risk groups. The K-M curve indicated that the high-risk group had a lower OS compared with the low-rank group (P=0.015) ().

When evaluating the OS-related prediction, the 1-, 3-, and 5-year AUCs of the developed gene signature were 0.970, 0.914, and 0.866, respectively, all of which were significant (). show the distribution of OS-related prediction models for patients in TCGA dataset, and it was evident that the OS-related prediction model were able to successfully differentiate between high-risk and low-risk groups.

shows the correlations between the 4 genes ALOX12, AKR1C2, SRC, and BNIP3 in the high- and the low-risk groups. In order to explore the correlations between these 4 genes and the prognosis of PCa patients, survival analysis was performed. High expression of SRC was significantly correlated with poorer OS in the K-M curve (P=0.00358) (). In addition, the genes ALOX12 (P=0.07935), AKR1C2 (P=0.19687), and BNIP3 (P=0.19119) were not significantly correlated with the OS rate of PCa patients ().

Figure 4

FRGs associated with the OS of prostate cancer patients. FRGs, ferroptosis-related genes; OS, overall survival.

The univariate Cox regression method was first to use to identify the 249 FRGs that were associated with DFS (). Hub genes were then selected using the LASSO method to minimize the risk of overfitting (). The multivariate Cox analysis was performed to determine the real DFS-related factors and finally identified a prognostic panel of 8 genes (ATG16L1, AKR1C2, MIOX, ALOX12, DRD4, IDH1, RELA, CHMP6), and a prediction model was established: DFS-related prediction model = (0.32010082 × expression value of ATG16L1) + (0.75803321 × expression value of AKR1C2) + (0.76159814 × expression value of MIOX) + (2.77524618 × expression value of ALOX12) + (0.97252682 × ex1-9pression value of DRD4) + (0.07452585 × expression value of IDH1) + (1.74889181 × expression value of RELA) + (−1.10159123 × expression value of CHMP6). Subsequently, we divided the 500 PCa cases into high- and low-risk groups according to the median of the DFS-related prediction model. The K-M plot indicated that the DFS rate of patients in the high-risk group was significantly lower than that of patients in the low-risk group (P=0.01) ().

Figure 5

DFS-related prognostic model of prostate cancer patients. (A) Analysis of 16 ferroptosis-related genes associated with DFS in PCa patients by univariate Cox regression. (B,C) A 37 mRNA signature was constructed using LASSO Cox regression. (D) Kaplan-Meier plot indicates that the DFS rate of patients in the high-risk group is significantly lower than that of the patients in the low-risk group. (E) ROC subject curve: the curve of AUC at 1 year, 3 years and 5 years, where AUC >0.75 predicts the results with statistical significance. (F,G) DFS-related prediction model distribution of patients in TCGA dataset. (H) Heat map showing the correlation of 8 genes between the high- and low-risk groups. AUC, area under the ROC curve; DFS, disease-free survival; PCa, prostate cancer; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Moreover, the AUC for 1 year, 3 years, and 5 years DFS were 0.980, 0.954, and 0.990, respectively, which were meaningful classification results (). show the distribution of DFS-related prediction models for patients in TCGA dataset. shows the correlations of 8 genes between high-risk and low-risk groups.

In order to explore the correlations between these 8 genes and the DFS rate of PCa patients, they were subjected to survival analysis. The results showed that ATG16L1 (P=0.06234), AKR1C2 (P=0.32121), MIOX (P=0.86353), ALOX12 (P=0.11591), DRD4 (P=0.31457), IDH1 (P=0.55412), RELA (P=0.31046), and CHMP6 (P=0.09613) were not significantly associated with the DFS rate of PCa patients ().

Figure 6

FRGs associated with the DFS of prostate cancer patients. FRGs, ferroptosis-related genes; DFS, disease-free survival.

Validation of predictive models using clinicopathological parameters from high- and low-risk patients

The OS-related prediction model showed statistically significant higher values for T3–4 than for T1–2 (P=0.0057), while no differences in OS-related prediction model values were observed between N0 and N1 (P=0.283), age ≥65 and age <65 (P=0.94), and Gleason score ≥7 and <7 (P=0.756) ().

Figure 7

Clinicopathological significance of OS-related prognostic models in prostate cancer. OS, overall survival.

The DFS-related prediction model was significantly different between T3–4 than T1–2 (P<0.001), N0 than N1 stage (P=0.0136); while no differences in DFS-related prediction model values were inspected between gleason score ≥7 and <7 (P=0.32), age ≥65 and age <65 (P=0.986) ().

Figure 8

Clinicopathological significance of the DFS-related prognostic model in prostate cancer. DFS, disease-free survival.

Ferroptosis subtypes were highly associated with pathway enrichment and immune infiltration among PCa risk groups

We examined the specific correlation of biological behaviors between PCa patients in the high- and low-risk groups, with the help of KEGG dataset to perform ssGSEA through the R package GSVA. As shown in , the p53 signaling pathway being the most enriched pathway in the high-risk group, and endocytosis being the most enriched pathway in the low-risk group.

Figure 9

Pathway enrichment and immune infiltration of ferroptosis subtypes differ among prostate cancer risk groups. (A) Enrichment analysis of metabolic pathways between the high- and low-risk groups. (B) Heat map showing immune infiltration based on CIBERSORT algorithms. (C) Heat map showing immune infiltration based on ssGSEA algorithms. (D) Boxplots showing differences in immune cell infiltration between the 2 ferroptosis high- and low-risk groups. (E,F) Bar plot and boxplot showing the relative abundance of ferroptotic immune cell infiltration in 500 prostate cancer patients. (G) Violin plot showing the ratio of StromalScore, ImmuneScore, and ESTIMATEScore between high-risk and low-risk groups. **, P<0.01; ns, P>0.05. GSVA, Gene set variation analysis; ssGSEA, single-sample gene set enrichment analysis.

In addition, the CIBERSORT algorithm was used to visualize the relative abundance of 12 immune-infiltrating cell subsets in the dataset by generating a heat map (). We found significant differences in ICI between the high- and low-risk groups (). M2-type macrophages (P=0.007) and neutrophils (P=0.024) were significantly enriched in the high-risk group, and CD4-activated memory T cells were significantly enriched in the low-risk group (P=0.017). Next, the abundance of immune infiltration of 22 immune cell subtypes in PCa tissues was assessed by a bar plot and boxplot. It was found that the relative numbers of resting memory CD4 T cells and plasma cells were the largest, while naive CD4 T cells had the lowest ratio (). Using the Estimate algorithm, a violin plot was used to show the ratio of StromalScore, ImmuneScore, and ESTIMATEScore between the high-risk group and low-risk group, and there was an apparent difference in StromalScore between the different risk groups in the immune matrix of the tumor microenvironment (P=0.0021) ().

Discussion

Plenty of evidence suggests that FRGs play a role in anti-proliferation and invasion through multiple signaling pathways (13). Ferroptosis is an iron oxide-dependent selective regulator of cell death (RCD) characterized by accumulating iron-dependent extensive lipid peroxidation (14). Because of its distinct biological behavior from apoptosis, necrosis, and autophagy, ferroptosis can play a role in many biochemical processes such as anabolism and redox processes, which was also the biological process of the outcome of most malignant tumors (15). So far, there have been several studies on the role of FRGs in PCa. Liu et al. established a model with 7 FRGs and concluded that high-risk PCa patients were positively correlated with tumor mutational burden (TMB) and bicalutamide resistance (16). Lv et al. established a prediction model including 9 FRGs and calculated that a high-risk score was positively correlated with the Ki67 expression level, iron potential index, and BCR risk in PCa patients, which illustrating that FRGs are associated with prognosis in PCa patients (17). With the advancement of accurate follow-up and individualized treatment of tumors, cancer treatment has entered a new era based on genes, that is, replacing tumor types based on tumor genes or molecular profiles (18). It is crucial to avoid overdiagnosis and overtreatment of tumors (19). In this study, FRG models were established to predict OS and DFS among PCa patients, and the immune infiltration landscape was systematically analyzed in an effort to perform a transcriptomic molecular assessment.

Enrichment analysis of the differential expression of FRGs was performed. The 249 RNA-seq expression data were collected from the TCGA-PRAD and the clinical data of 500 patients. The GO results showed that the most enriched terms were certain metabolic pathways, which is consistent with the pathophysiological mechanism of ferroptosis confirmed in previous studies (20). KEGG results with higher enrichment included microRNAs in cancer, ferroptosis, arachidonic acid metabolism, and thyroid hormone synthesis. Interferon γ (IFNγ) binding to arachidonic acid (AA) induced immunogenic tumor ferroptosis as a mode of action in CD8+ T cell (CTL)-mediated tumor killing (21). Study by Mishima et al. suggested that reverse T3, a physiologically inactive form of triiodothyronine (T3), also prevents ferroptosis like other thyroid hormones (22).

Based on the difference analysis and univariate and LASSO regression analysis results of tumor and adjacent tissues, 31 FRGs associated with OS and 16 FRGs associated with DFS were obtained. Finally, the prediction models of OS and DFS were established by multivariate Cox regression analysis, which proved to be independent prognostic indexes for PCa patients.

Among ALOX12, AKR1C2, SRC, and BNIP3, and the 4 genes involved in the establishment of the OS-related prognostic model, high expression of SRC (P=0.04) was significantly correlated with poor OS in the K-M curve. ALOX12 plays a role in tumorigenesis by promoting cellular proliferation, amino acid metabolism and lipid synthesis to trigger inflammatory responses (23). the expression of ALOX12 is significantly downregulated and inhibits the proliferation of pancreatic adenocarcinoma (PAAD) cells (24), as well as related to the stage and grade of PCa (25). In prostate cells, AKR1C2 acts as a 3-ketosteroid reductase to reduce dihydrotestosterone (DHT) levels and block the activation of androgen receptor (26). Selective deletion of AKR1C2 in PCa may promote clonal expansion of tumor cells by reducing DHT metabolism and enhancing androgen-dependent cell proliferation (27). SRC is a proto-oncogenic non-receptor tyrosine kinase that plays a role in tumor cell proliferation, survival, and invasion. SRC is one of the most active kinases in phosphoproteomics of aggressive PCa (28). Chakraborty et al. demonstrated that the proto-oncogene SRC reduced primary tumor growth and metastasis by reprogramming glucose metabolism through the encoded non-receptor tyrosine kinase c-Src (29). SRC activation was found in 28% of men who developed castration-resistant PCa after androgen deprivation therapy (30). BNIP3, regarding as a member of the BH3 family, locate mitochondria and is associated with apoptosis, autophagy, or mitophagy. Elevated expression of BNIP3 in human PCa predicts poorer prognosis (31).

Among the 8 genes involved in establishing the DFS prognostic model, namely ATG16L1, AKR1C2, MIOX, ALOX12, DRD4, IDH1, RELA, and CHMP6, Hurley et al. confirmed that charged multivesicular body protein 6 (CHMP6) is ingredient of endosomal sorting complexes required for transport (ESCRT), which is required for many cellular processes involving membrane remodeling (32). Gene knockdown of CHMP6 significantly enhances tau aggregation, resulting in damage to the endolysosomal membrane, consistent with a role for the ESCRT pathway in endolysosomal membrane repair (33). In addition, mutations in CHMP6 components have been shown to be associated with degeneration of the nervous system (34). However, so far, there is no research on CHMP6 in the field of cancer. Our study shows that its high expression is significantly correlated with better DFS in PCa, which needed further exploration. DRD4 is one of the subtypes of the dopamine receptor gene in G protein-coupled receptors (GPCRs) (35), containing a variable number of 48 bp tandem repeats (VNTRs) in the third exon, and its homozygous and heterozygous forms are associated with better and worse prognosis in PCa, respectively (36). Its expression level also correlated with prostate specific antigen (PSA) levels in PCa (37). Related studies have shown that different mutation types of the metabolic enzyme isocitrate dehydrogenase 1 (IDH1) induce ferroptosis through different pathways. Pathways that inhibit ferroptosis include mutant IDH1 to produce the oncogenic metabolite D-2-hydroxyglutarate (D-2-HG) and deletion of mutant IDH1 alleles in IDH1 heterozygous tumor cells. The pathway that promotes the accumulation of lipid ROS, which in turn promotes ferroptosis, has ectopic expression of mutant IDH1 (38). High ATG16L1 expression is a poor prognostic marker in multiple cancer types, and its knockout resensitizes lung cancer cells to cancer therapy after TGF-β-induced resistance (39). However, in a clinical study involving 458 PCa patients, the results showed that high expression of ATG16L1 was associated with lower tumor aggressiveness and better prognosis (40). Therefore, it is necessary to further explore and study the correlation between ATG16L1 expression and the prognosis of PCa patients.

The prediction model was verified by clinical data such as age, Gleason score, T stage, and N stage of PCa patients. The results showed that the OS- and DFS-prediction models were correlated with the T stage score, which was higher in T3–4 than in T1–2 (P=0.0057, P=0.000697). Moreover, the N0 stage was higher than the N1 stage in the DFS prediction model (P=0.0136), with a significant difference. T stage is an important prognostic indicator for PCa patients and is positively correlated with BCR (41). As for PCa, especially in the stage of CRPC, the occurrence of lymph node metastasis is an important indicator for judging the prognosis of patients (42).

Special biological properties of the tumor cell transcriptome, such as metabolic reorganization and immune evasion, make it resistant to multiple chemotherapeutic drugs (43). This not only affects the treatment effect, but also jeopardizes the prognosis of patients (44). We evaluated the immune infiltration landscape of FRGs in PCa and attempted to comprehensively analyze the relationship between the degree of ICI and the activity of immune-related functions in the high- and low-risk groups. Through GSEA, we found that P53 gene mutations and the pentose phosphate pathway were enriched in high-risk groups of PCa. TP53 is one of the most frequently mutated suppressor genes in PCa, and its inactivation is a mid-to-late event in PCa progression, mainly associated with mCRPC (45). Strohecker et al. described a putative stimulator of autophagy, 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 4 (PFKFB4), that drives flux through the pentose phosphate pathway. Knockdown of PFKFB4 in PCa cells increases p62 and ROS, but surprisingly increases autophagic flux (46). Interestingly, immune genes were equally present in the high- and low-risk groups, regardless of their enrichment status, and these results are helpful for establishing prognostic risk assessment models.

In the risk score model derived from the CIBERSORT algorithm, we found that the relative numbers of resting memory CD4+ T cells and plasma cells were the largest, while the proportion of naive CD4+ T cells was the lowest, which is highly consistent with current research on ICI. Wang et al. reported that UBASH3B (STS1) is an important gene that negatively regulates T cell receptor signaling in activated T lymphocytes involved in cancer and is associated with resting CD4+ memory T cells or de novo PCa therapeutic targets (47). IL-6 produced by prostate infection with trichomonas vaginalis impacts on the tumor microenvironment by promoting PCa cell progression after M2 macrophage polarization induction (48). The number of activated CD4+ memory T cells in the high TMB group was significantly higher than that in the low TMB group, and the TMB was associated with BCR-free survival and the infiltration of activated CD4+ memory T cells in the immune microenvironment (49). The above is consistent with our results showing that M2-type macrophages were significantly enriched in the high-risk group, while CD4-activated memory T cells were highly expressed in the low-risk group. Taken together, the risk score obtained in this study is associated with the immunosuppressive microenvironment of the tumor. Finally, through the Estimate algorithm, it could be concluded that the immune matrix score of the high-risk group was higher than that of the low-risk group. PCa is considered an immunogenic tumor. A phase II clinical trial jointly conducted by the Institute of Cancer Research London and Royal Marsden Foundation Trust suggests that some patients with advanced PCa do respond, especially those with mutations in DNA repair genes that respond well to immunotherapy (50).

In this study, we established FRG-related OS and DFS-related prediction models of PCa, and for the first time proposed that the expression of CHMP6 was associated with prognosis. It may be related to the mechanism of endosomal sorting and endolysosomal membrane repair. As for ATG16L1, it has been reported that its high expression is associated with poor cancer prognosis, but previous studies have clarified that this gene can be used as an indicator of better prognosis in PCa. Therefore, analyzing and determining the clinical significance of the conduction pathway of CHMP6 and ATG16L1 in PCa is necessary as it may become a new research and therapeutic target. However, there are many shortcomings in this study. Firstly, this study is a retrospective data analysis. Without prospective clinical trials, strong evidence cannot be provided. Secondly, the data are obtained from public platform resources, such as TCGA and KEGG databases. Since no external validation was done, the accuracy was limited. There was no data from our center, and the clinicopathological data of patients in this clinical research center can be added in subsequent studies. Undoubtedly, there is still a need for a large number of validation studies in the future.

The article’s supplementary files as