The transcriptome data (TCGA-PRAD) on gene expression matrix and relapse information of PCa were extracted from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). The TCGA-PRAD dataset (access time: November 20th, 2024) included 502 PCa tissue samples (397 samples with relapse information) and 52 paracancer (control) tissue samples. Among the 397 samples, 70% samples (278 samples) were employed as the training set (TCGA-PRAD-train), and 30% samples (119 samples) were employed as the validation set (TCGA-PRAD-validation). The single-cell RNA sequencing (scRNA-seq) data (GSE141445) of PCa was scoured from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE141445 (GPL24676, access time: November 20th, 2024) included 13 PCa tissue samples. Validation dataset GSE116918 contained 248 samples with complete recurrence information. The MR data (EBI-A-GCST90018905) of PCa and eQTL data were acquired from the Integrative Epidemiology Unit (IEU) Open Genome-wide Association Study (GWAS) database (https://gwas.mrcieu.ac.uk/). The “EBI-A-GCST90018905” included 24,119,306 SNPs from 211,227 Europeans (case: control=11,599: 199,628). The 423 MCD-related genes (MCDRGs) used in the study were obtained by merging and removing duplicates from three myeloid cell differentiation-related datasets: GOBP_NEGATIVE_REGULATION_OF_MYELOID_CELL_DIFFERENTIATION (95 genes), GOBP_POSITIVE_REGULATION_OF_MYELOID_CELL_DIFFERENTIATION (104 genes), and GOBP_MYELOID_CELL_DIFFERENTIATION (423 genes), which were downloaded from the MCD related genes (MCDRGs) were acquired from Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/human/search.jsp) and relevant references (20) ( Supplementary Table 1 ).

To acquire differentially expressed genes (DEGs) between PCa and control samples in TCGA-PRAD, “DESeq2” package (v 3.4.1) was carried out (PCa vs control) (|log₂Fold Change (FC)| > 0.5, P < 0.05) (21). On the basis of the log₂FC value, DEGs were visualized and the top 10 up/down-regulated gene names were labeled by the volcano plot utilizing “ggplot2” package (v 3.4.1) (22). Similarly, the expressions of the top 10 up/down-regulated genes between PCa and control groups were displayed by the heat plot utilizing “ComplexHeatmap” package (v 2.14.0) (23).

DEGs were intersected with MCDRGs to acquire candidate genes via “ggvenn” package (v 1.7.3) (24). To probe the organic activities and signal pathways involved in candidate genes, based on background gene set of “org.Hs.eg.db” package (v 3.16.0) (25), Gene Ontology (GO) (P < 0.05) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (P < 0.05) enrichment analyses were carried out utilizing “clusterProfiler” package (v 4.7.1.3) (26). Subsequently, to determine the interactions of proteins encoded by candidate genes, candidate genes were uploaded to STRING database (https://cn.string-db.org/) (interaction score > 0.4). The results were imported into Cytoscape software (v 3.9.1) and Protein-Protein Interaction (PPI) network was constructed (27).

To probe the causal dependence between candidate genes and PCa and acquire candidate prognostic genes, MR analysis was performed with PCa as outcome event and candidate genes as exposure factors utilizing “TwoSampleMR” package (v 0.6.1) (28). The MR consisted of 3 assumptions: (1) instrumental variables (IVs) were linked with exposure factors; (2) IVs could only influence outcomes by exposure factors; (3) IVs were not linked with potential confounders. To obtain effective IVs, exposure factors and outcome event were read and IVs were screened via “extract_instruments” function. The screening criteria were as follows: (1) IVs strikingly linked with exposure factors (P < 5×10⁻⁶); (2) IVs that exhibited linkage disequilibrium were removed with R² = 0.001, kb=10, and clump=TRUE; (3) IVs that were strikingly linked with outcome were removed with rsq=0.8 and proxies=TRUE; (4) IVs whose F statistic < 10 were removed ( F = ( N − K − 1 ) K / R 2 1 − R 2 , R² represented the cumulative explanatory variance of SNPs, N represented the number of the samples, and K represented the number of SNPs). Then, the effect alleles and effect sizes were unified via “harmonCIe_data” function. The 5 algorithms of the “mr” function were utilized to conduct MR analysis for each exposure factor and outcome, which included Inverse Variance Weighted (IVW) (29), Weighted Mode (30), MR Egger (31), Simple Mode (28), and Weighted Median (32). MR analysis mainly relied on IVW results. Exposure factors with SNP > 2 and P < 0.05 were considered as the exposure factors that had a causal relationship with PCa. The odds ratio (OR) > 1 suggested that the exposure factor was a contributing factor to the risk of PCa, while an OR < 1 indicated that it was a protective factor. Notably, the scatter plot was drawn to further identify the correlation between exposure factors and outcome in combination with SNP-exposure effects and SNP-outcome effects via “mr_scatter_plot” function. The forest plot was drawn to evaluate the diagnostic power of the estimated exposure factors of each SNP site on the outcome via “mr_forest_plot” function. The funnel plot was drawn to judge whether the analysis was random and adhered to Mendel’s second law random grouping in combination with the β and standard error (SE) of each IV via “mr_funnel_plot” function.

Moreover, to evaluate the reliability of MR analysis results, sensitivity analysis was applied. Heterogeneity test was performed via “mr_heterogeneity” function (P > 0.05, Cochran’s Q test). Horizontal pleiotropy test was performed to find out whether confounding factors existed via “mr_pleiotropy_test” and “MR-Egger” functions (P > 0.05). To observe whether the SNPs of each IV caused considerable changes in the outcome, Leave-one-out (LOO) analysis was performed via “mr_leaveoneout” function.

To verify that the results of the forward analysis were not interfered by reverse causality and determine the validity of the causal sequence between the outcome and the exposure factors, Steiger test was carried out via “steiger_filtering” function (Steiger-dir=TRUE, P < 0.05).

At last, the candidate genes examined by sensitivity analysis and Steiger test were specified as candidate prognostic genes.

To identify prognostic genes, in TCGA-PRAD-train, univariate Cox regression analysis was performed on the basis of candidate prognostic genes via “survival” package (v 3.5-3) (33)(P < 0.2) (34, 35) and the result was presented via “forestplot” package (v 2.0.1) (36). Genes with consistent hazard ratios (HRs) or ORs in both univariate Cox regression analysis and MR analysis were utilized for proportional hazards (PH) assumption test via “cox.zph” function (P > 0.05). After that, Least Absolute Shrinkage and Selection Operator (LASSO) was implemented utilizing “glmnet” package (v 4.1.4) (37)(10-fold cross-validation). Finally, genes that had passed through the above analyses in sequence were defined as prognostic genes to construct a risk model.

In TCGA-PRAD-train, on the basis of the expression of prognostic genes and the risk coefficients gained from LASSO regression, PCa patients were scored with the following formula.

Expr represented the expression level of each prognostic gene and coef signified the risk coefficient of each prognostic gene. Notably, according to the median of risk score, PCa patients were categorized into high risk group (HRG) and low risk group (LRG) and risk score distribution and survival status of patients were displayed. According to relapse time and relapse status of PCa patients, the “survminer” package (v 0.4.9) (38) was implemented to draw and compare survival curves of HRG and LRG (P < 0.05). The Area Under Curve (AUC) values of receiver operating characteristic (ROC) curves at 1, 2, and 3 years were employed to evaluate the accuracy of risk model utilizing “survivalROC” package (v 1.18.0) (39) (AUC > 0.6). The validation was performed in the GSE116918 dataset. Additionally, the heat plot was drawn to display the expression levels of prognostic genes in HRG and LRG. Using the same method, 2 risk models were constructed in TCGA-PRAD-validation and all samples with relapse information of TCGA-PRAD to validate the accuracy of the above-mentioned model.

To evaluate the ability of prognostic genes to predict PCa recurrence rates, the “regplot” package (v 1.1) (40) was employed to build the nomogram model for PCa in the TCGA-PRAD-train. Each prognostic gene was scored separately, and the scores were added together to obtain the total scores. The higher the total scores, the higher the recurrence rate of the patient. Calibration curves built via “rms” package (v 6.5.0) (41) and ROC curves built via “survivalROC” package (v 1.18.0) (AUC > 0.6) were applied to compute the effectiveness of the nomogram model in clinical prediction at 1, 2, and 3 years for PCa.

To estimate the tumor purity and the status of stromal and immune cells in the malignant tumor tissues within the tumor microenvironment of PCa patients. The data were normalized using the R package estimate (v 1.0.13, https://R-Forge.R-project.org/projects/estimate/), and the stromal score, immune score, and ESTIMATE score were calculated using the ESTIMATE algorithm. Wilcoxon test was employed to contrast the scores mentioned above between HRG and LRG (P < 0.05). For further evaluation of the situation of immune cells in the development process of PCa, based on single sample Gene Set Enrichment Analysis (ssGSEA) algorithm, in the TCGA-PRAD-train, “GSVA” package (v 1.46.0) (42) was applied to compute the enrichment scores of 28 immune cells (43) in the HRG and LRG. Wilcoxon test was implemented to contrast differences in immune cell infiltration between HRG and LRG (P < 0.05). What’s more, to understand the link between prognostic genes and differential immune cells, in the TCGA-PRAD-train, “psych” package (v 2.2.9) (44) was performed (|correlation (cor)| > 0.3, P < 0.05).

To determine the biological pathways involved in the occurrence and development of PCa in the HRG and LRG, in the TCGA-PRAD-train, “DESeq2” package (v 3.4.1) was employed to perform differential expression analysis between HRG and LRG and log₂FC values were calculated. The log₂FC values were sorted in descending order. Then, Gene Set Enrichment Analysis (GSEA) was performed via “clusterProfiler” package (v 4.7.1.3) (|Normalized Enrichment Score (NES)| > 1, P < 0.05). The top 5 pathways with notable P-values were presented. Besides, GeneMANIA database (http://genemania.org) was applied to predict the genes related to the functions of prognostic genes and their involved functions.

To probe the drug sensitivity of the HRG and LRG, drugs related to tumors were obtained from Genomics of Drug Sensitivity in Cancer 2 (GDSC2) database (https://www.cancerrxgene.org/). In the TCGA-PRAD-train, “pRRophetic” package (v 0.5) (45) was employed to determine the half maximal inhibitory concentration (IC50) of each tumor sample. Wilcoxon test was employed to contrast the differences in drug sensitivity between the HRG and LRG (P < 0.05). The result was displayed via “ggplot2” package (v 3.4.1).

To explore the upstream regulatory factors of prognostic genes and their interaction relationships, based on NetworkAnalyst online website (https://www.networkanalyst.ca/NetworkAnalyst/uploads/ListUploadView.xhtml), Transcription Factors (TFs) were predicted by TRRUST database (https://www.grnpedia.org/trrust/) and miRNAs were predicted by miRTarBase database (v 9.0) (https://mirtarbase.cuhk.edu.cn/) and TarBase (v 9.0) (https://dianalab.e-ce.uth.gr/tarbasev9). Finally, the TF-mRNA-miRNA network was constructed utilizing Cytoscape software (v 3.9.1).

All scRNA-seq analyses were performed via “Seurat” package (v 5.0.1) (46). To ensure the accuracy and reliability of scRNA-seq data, in the GSE141445, the data was filtered via “PercentageFeatureSet” function. The screening criteria were as follows: (1) the number of genes in the cells ranged from 200 to 3,000; (2) genes with an expression level ranging from 200 to 4,531 and covered by at least 3 cells were retained. Then, the data after filtering were normalized via “NormalizeData” function. The 2,000 highly variable genes (HVGs) were acquired and the top 10 most mutated genes were labeled via “FindVariableFeatures” function and the result was presented via “LablePoints” function. Moreover, the samples of GSE141445 were subjected to normalization processing via “Scale Data” function. The HVGs were subjected to principal component analysis (PCA) via “runPCA” function. The P-value when the principal components (PCs) ranged from 1 to 30 was calculated via “Jackstraw” function. The values of the sudden variance drops when different values were taken for the PCs were calculated and the result was displayed via “Elbowplot” function. When P < 0.05, the PCs at the inflection point in the variance elbow plot were used for subsequent analysis. Additionally, unsupervised clustering analysis was applied via “FindClusters” and “FindNeighbors” functions (resolution=0.1). Notably, cells were clustered and result was visualized via “RunTSNE” function. Furthermore, to determine the cell types of the cell clusters, cell clusters were commented on the basis of marker genes (47). The expression patterns of marker genes in different cell clusters and cells were displayed through bubble plots.

To gain key cells associated with PCa development, in the GSE141445, the manifestation of prognostic genes in cells were displayed through bubble plots. The cells with the highest gene expression were labeled as key cells. Besides, the distribution of prognostic genes in cells was also presented.

In the GSE141445, the “CellChat” package (v 1.6.1) (48) was employed to explore the communication networks among different cells. The pairing of ligands and receptors among cells was presented through bubble diagrams using “ggplot2” (v 3.4.1). Furthermore, to investigate the differentiation of key cells, key cells were subjected to secondary clustering using the same method as that in 2.11. Then, the developmental trajectories of key cells were analyzed via the “Monocle” package (v 2.22.0) (49). The “DDRTree” package (v 0.1.5) (50) was employed to draw cell trajectory diagram. Finally, the expression of prognostic genes during the differentiation process of key cells was also explored.

To clarify the manifestation of prognostic genes, in the PCa tissue and paracancer tissue samples of TCGA-PRAD, Wilcoxon test was employed to compare the differences in manifestation of prognostic genes between PCa samples and control samples (P < 0.05). Subsequently, to verify the accuracy of the above results, RT-qPCR was performed. A total of 5 pairs of tissue samples were obtained at Tangdu hospital, including 5 PCa and 5 paracancer (control). Informed consent was obtained from all patients for the use of PCa tissue samples in this study. The informed consent form needed to be signed and filled out by all participants, while the ethical approval agency was The Clinical Ethics Committee of Tangdu Hospital of Air Force Medical University (permission number: TDLL-KY-202405-18). Then, total RNA of 5 pairs of tissue samples was extracted by TRizol reagent (Vazyme, Nanjing, Jiangsu). The RNA concentrations were computed by NanoPhotometer N50. Subsequently, mRNA was converted to cDNA by Hifair^® III 1st Strand cDNA Synthesis SuperMix for qPCR test kit (Yeasen, Shanghai). Finally, RT-qPCR was carried out. The primers of reaction reagents, reaction conditions and genes were arranged in Supplementary Table 2 . The internal reference gene was GAPDH, which was employed to normalize the results. The expression levels of prognostic genes were calculated by 2^-ΔΔCt. The results were calculated by GraphPad Prism software (v 5.0) (51).

Cells were lysed in RIPA buffer (P0013B, Beyotime) to extract proteins, and protein concentrations were determined using a BCA assay kit (P0012, Beyotime). Samples were separated by SDS-PAGE (12%, Invitrogen) and transferred to PVDF membranes, which were then blocked with BSA for 2.5 h. After washing, membranes were incubated overnight at 4°C with primary antibodies, including anti-BMP2 (ab284387, 1:1000, Abcam) and anti-FASN (ab128870, 1:1000, Abcam). The next day, membranes were washed three times and incubated with HRP-conjugated goat anti-rabbit IgG secondary antibody (G-21234, 1:5000, Invitrogen) at room temperature for 2 h. Chemiluminescent substrate (ECL) was applied evenly to the membranes, and signals were captured using a ChemiDoc XRS+ imaging system (Bio-Rad, USA). Band intensities were analyzed with ImageJ, and protein levels were quantified as the ratio of the target band intensity to that of β-actin (ab272085, 1:1000, Abcam).

Vector, OE-BMP2, sh-NC, and sh-FASN cells were seeded and allowed to adhere for 24 h. Cells were washed twice with phosphate-buffered saline (PBS) and fixed with 4% paraformaldehyde for 20 min. Permeabilization was performed with 0.3% Triton X-100 for 15 min, followed by blocking with 6% donkey serum for 30 min. Cells were incubated overnight at 4°C with primary antibody against Ki67 (1:200 in PBS). The next day, cells were incubated with a fluorescent secondary antibody (1:500 in PBS) for 1 h at 37°C in the dark. After three PBS washes, nuclei were counterstained with DAPI for 10 min. Ki67 expression was assessed using a fluorescence microscope.

Vector, OE-BMP2, sh-NC, and sh-FASN cells (1×10^6 cells/well) were seeded in 6-well plates and allowed to adhere for 24 h. Once confluence exceeded 90%, a straight scratch perpendicular to the horizontal axis was made using a sterile 20-µL pipette tip. To remove debris, cells were rinsed with serum-free DMEM. Plates were incubated, and scratch closure was monitored at 24 h and 48 h after wounding under a high-power microscope. Migrating cells within the scratch area were quantified using ImageJ.

Human prostate cancer PC3 cells (ATCC^® CRL-1435™) were maintained in RPMI-1640 medium (Gibco) supplemented with 10% fetal bovine serum (FBS; Corning) and 1% penicillin/streptomycin (HyClone) at 37°C with 5% CO₂. Lentiviral particles expressing shRNA targeting human fatty acid synthase (FASN; target sequence: shRNA1:5′-GCATGGAGCGTATCTGTGAGAACTCGAGTTCTCACAGATACGCTCCATGTTTTTT-3′ shRNA2: 5′-GCTACGACTACGGCCCTCATTCTCGAGAATGAGGGCCGTAGTCGTAGCTTTTTT-3′) or non-targeting scramble control (shRNA-NC: 5′-GCACCCAGTCCGCCCTGAGCAAATTCAAGAGATTTGCTCAGGGCGGACTGGGTGCTTTTT-3′) in pLKO.1 vector were packaged in HEK293T cells using psPAX2 and pMD2.G plasmids. Viral supernatants were concentrated by ultracentrifugation (26,000 × g, 2.5 h) and titers determined via qPCR.

The human BMP2 ORF (NM_001200) was cloned into pLV-CMV-3flag-zsgreen vector via NotI and NheI site. Lentivirus production and titration followed procedures (titer ≥5×10⁷ TU/mL). PC3 cells were transduced at MOI=20 with polybrene supplementation.

Cells were seeded in 6-well plates (5×105 cells/well) and incubated for 24 h to reach 50–60% confluency. Transduction was performed by replacing medium with viral suspension (MOI=20) containing 8 μg/mL polybrene (Sigma).

R software (v 3.9.1) was implemented to apply bioinformatics analyses. Wilcoxon test and t test were utilized to compare the disparities between 2 groups (P < 0.05).

By differential expression analysis, 13,091 DEGs were acquired, comprising 7,935 DEGs with up-regulated expression and 5,156 DEGs with down-regulated expression. All DEGs and the top 10 up/down-regulated gene names were presented ( Figure 1A ). Besides, the expressions of the top 10 up/down-regulated genes in PCa and control groups were also displayed ( Figure 1B ). After that, 141 candidate genes were acquired ( Figure 1C ). GO analysis indicated that candidate genes were enriched in 1,758 functions, including 1,619 biological process (BP) items such as MCD, 37 cellular component (CC) items such as chromosome centromeric core domain, and 102 molecular function (MF) items such as cytokine receptor binding ( Figure 1D , Supplementary Table 3 ). KEGG analysis suggested that candidate genes were enriched in 85 pathways such as osteoclast differentiation ( Figure 1E , Supplementary Table 4 ). Subsequently, PPI network indicated that there were interactions among proteins encoded by 122 candidate genes such as MYC, H4C11, and MMP9 ( Figure 1F ).

In order to investigate the causal relationship between DE-MCDRGs and PCa outcomes, MR analysis was conducted using 141 candidate genes as exposure factors and PCa as the outcome. Firstly, the 46 candidate genes had a substantially causal relationship with PCa based on the IVW method (P < 0.05), including 23 risk factors (OR > 1) and 23 protective factors (OR < 1) ( Table 1 ). The positive slope of the scatter plot indicated that the gene was a risk factor, while a negative slope indicated that it was a protective factor. The intercept close to 0 suggested the absence of confounding factors ( Supplementary Figure 1 ). For protective factors, the effect size of each SNP was less than 0. In contrast, for risk factors, the effect size of each SNP was greater than 0 ( Supplementary Figure 2 ). The symmetrical arrangement of SNPs around each exposure indicated that MR affirmed Mendel’s second law of randomization ( Supplementary Figure 3 ). So 46 candidate genes were used for subsequent analysis. Notably, among the 46 candidate genes, heterogeneity was calculated for 35 of them. The P-values of 34 candidate genes were greater than 0.05, indicating the absence of heterogeneity. Since the P-value of SLC4A1 was less than 0.05, SLC4A1 was analyzed with random effects IVW ( Supplementary Table 5 ). Meanwhile, among the 46 candidate genes, the P-values of 33 candidate genes were greater than 0.05, suggesting the lack of horizontal pleiotropy and the reliability of the results. However, the P-values of the remaining 13 candidate genes were less than 0.05, suggesting the presence of horizontal pleiotropy. To ensure the reliability of the results, these 13 candidate genes were excluded ( Supplementary Table 6 ). LOO analysis demonstrated that the removal of any SNP had minimal impact on the results ( Supplementary Figure 4 ). Steiger test indicated that only the P-values of 33 candidate genes were all less than 0.05 and the Steiger-dir values were all TRUE ( Supplementary Table 7 ). Finally, taking all of the above screening results into account, only 23 candidate genes could be labeled as candidate prognostic genes.

In order to further screen for genes related to prognosis among the 23 candidate prognostic genes, univariate Cox regression analysis was conducted. The 9 candidate prognostic genes associated with PCa recurrence were gained utilizing univariate Cox regression analysis (P < 0.2), including 4 risk factors (HR > 1) and 5 protective factors (HR < 1) ( Figure 2A ). Notably, compared with the results of the MR analysis, only the OR/HR of 6 candidate prognostic genes was consistent. Besides, these 6 candidate prognostic genes were also identified by PH assumption test (P > 0.05) ( Figure 2B ). Among them, 5 candidate prognostic genes were identified by LASSO analysis (lambda.min=0.005517655) ( Figure 2C ). So NR3C1, BMP2, RACGAP1, TLR3, and FASN were considered as prognostic genes to build the risk model in TCGA-PRAD-train. On the basis of the median of the risk scores (1.111923), PCa patients were divided into a HRG (139 PCa patients) and a LRG (139 PCa patients). Patients’ risk score distribution ( Figure 2D ) and survival status ( Figure 2E ) were displayed. The survival curves indicated that the survival rate of the HRG was substantially lower than that of the LRG (P=0.00077) ( Figure 2F ). The AUC values of the ROC curves were all greater than 0.7, indicating that the risk model had good performance ( Figure 2G ). In addition, BMP2, RACGAP1, and FASN were prominently expressed in the HRG, while TLR3 and NR3C1 were prominently expressed in the LRG ( Figure 2H ).

In order to assess the accuracy of risk prediction, the ROC curve was used to calculate the AUC value based on the TCGA-PRAD training set and the GSE116918 dataset to evaluate the model’s effectiveness. In the TCGA-PRAD-validation, on the basis of the median of the risk scores (1.056659), PCa patients were separated into a HRG (59 PCa patients) and a LRG (60 PCa patients). Patients’ risk score distribution ( Figure 3A ) and survival status ( Figure 3B ) were also displayed. The survival curves indicated that HRG had a lower survival rate than LRG (P=0.033) ( Figure 3C ). As the AUC values were all above 0.7, it demonstrated that the risk model had excellent performance ( Figure 3D ). BMP2, RACGAP1, and FASN exhibited high expression in the HRG, with TLR3 and NR3C1 showing high expression in the LRG ( Figure 3E ). Similarly, in all samples with relapse information of TCGA-PRAD, on the basis of the median of the risk scores (1.085715), PCa patients were categorized into a HRG (198 PCa patients) and a LRG (199 PCa patients). The risk score distribution ( Figure 3F ) among patients and their survival status ( Figure 3G ) were likewise shown. According to the survival curves, the survival rate in the HRG was markedly lower than that in the LRG (P < 0.0001) ( Figure 3H ). Given that the AUC values all exceeded 0.7, it was evident that the risk model performed well ( Figure 3I ). High expression of BMP2, RACGAP1, and FASN was observed in the HRG, while high expression of TLR3 and NR3C1 was found in the LRG ( Figure 3J ). Based on the results of 3.4, the risk model was accurate. In the GSE116918 dataset, there was a significant survival difference between the high and low-risk groups (p=0.014), but the ROC curves for recurrence at 1, 2, and 3 years showed that the AUC values were all less than 0.6. This suggests that although the risk score was significantly associated with survival differences, its accuracy in predicting recurrence remains limited. Further validation and optimization may be needed, possibly by integrating other clinical indicators or more comprehensive models ( Supplementary Figure 5 ).

The nomogram indicated that prognostic genes could predict the return rate of PCa patients quite well, and the recurrence rate increased as the total score of the nomogram rose ( Figure 4A ). The slope of all calibration curves was close to 1, indicating the accuracy of the nomogram model ( Figure 4B ). The AUC values of ROC curves all exceeded 0.7, which further verified the accuracy of the nomogram model ( Figure 4C ). In conclusion, prognostic genes could be utilized to predict the recurrence rate of PCa patients.

In order to assess tumor purity, immune infiltration, and stromal infiltration in the malignant tumor tissue within the patient’s tumor microenvironment, immune infiltration analysis was conducted. Within the microenvironment specific to PCa tumors, immune scores, ESTIMATE scores, and the stromal scores in the LRG were all strikingly higher than those in the HRG (P < 0.05) ( Figure 5A ). Subsequently, the infiltration levels of 28 immune cells in the HRG and LRG were analyzed ( Figure 5B ). There were 18 immune cells that showed notable distinctions between the HRG and the LRG, and they were defined as differentially expressed immune cells (P < 0.05) ( Figure 5C , Supplementary Table 8 ). Except for the activated CD4 T cells, the extents of infiltration of the remaining cells in the LRG were all remarkably higher than those in the HRG (P < 0.05). TLR3 had the largest notable positive link with natural killer cells (cor=0.66, P < 0.0001) and FASN had the largest notable negative correlation with natural killer (NK) T cells (cor=-0.38, P < 0.0001) ( Figure 5D , Supplementary Table 9 ). In summary, the occurrence and development of PCa might be related to changes in the tumor microenvironment or some immune cells.

In order to determine the biological pathways involved in the development of PCa between the high and low-risk groups, GSEA enrichment analysis was conducted. By GSEA analysis, the HRG and LRG were strikingly enriched 41 pathways such as hcm, porphyrin and chlorophyll metabolism, pentose and glucuronate interconversions, and cell cycle ( Figure 6A , Supplementary Table 10 ). These pathways suggested that the risk of PCa might be related to some biological role. In order to obtain molecular targeted drugs corresponding to the genes in the high and low-risk groups, drug sensitivity analysis was conducted. Among 198 drugs, the HRG and LRG suggested notable disparities in sensitivity to 86 drugs (P < 0.05) ( Figure 6B , Supplementary Table 11 ). The IC₅₀ value of AZD8186 in the HRG was strikingly higher than that in the LRG, indicating that the LRG was more sensitive to this drug. On the contrary, the IC₅₀ value of ML323 in the HRG was remarkably lower than that in the LRG, indicating that the HRG was more responsive to this drug.

In order to identify other genes related to the function of prognostic genes, the GeneMANIA database was used to predict genes associated with the function of prognostic genes and the functions they are involved in. By GeneMANIA analysis, 20 genes related to the functions of prognostic genes were acquired such as BMPR1A and these genes were related to 7 functions such as response to BMP ( Figure 7A ). Furthermore, in order to explore the upstream regulatory factors and their interactions for the prognostic genes, a TF-mRNA-miRNA regulatory network was constructed. Two TFs and 32 miRNAs were predicted for FASN; 5 TFs and 32 miRNAs were predicted for NR3C1; 2 TFs and 32 miRNAs were predicted for RACGAP1; 4 TFs and 10 miRNAs were predicted for BMP2; 6 TFs and 10 miRNAs were predicted for TLR3 ( Figure 7B ). Then, the network was constructed to demonstrate the complex relationships between prognostic genes and regulatory molecules.

In order to ensure the accuracy, reliability, and interpretability of the single-cell data, quality control was performed on all scRNA-seq data. There were 24,391 genes and 36,424 cells before quality control ( Figure 8A ), and 24,391 genes and 34,571 cells after quality control ( Figure 8B ). Then, the 2,000 HVGs and the top 10 most varied gene names were displayed ( Figure 8C ). After dimensionality reduction, the first 30 PCs were used for clustering analysis ( Figure 8D ). Finally, the 14 cell clusters were acquired ( Figure 8E ). After annotating cell clusters with markers genes, there were 7 cells acquired, which included mast cells (TPSB2, MS4A2, TPSAB1), stroma cells (COL1A2, TAGLN, ACTA2), endothelial cells (PECAM1, VWF, ACKR1), T cells (GZMA, CD3E, CD3D), epithelial cells (KRT18, EPCAM, KRT19), myeloid cells (CD14, FCGR3A, CD163), and B cells (MS4A1, CD79A, DERL3) ( Figures 8F-H ). Among all the prognostic genes, FASN was prominently expressed in epithelial cells, so epithelial cells were considered as key cells ( Figure 8I ). NR3C1 was more abundantly distributed in T cells, and FASN was more abundantly distributed in epithelial cells ( Figure 8J ).

In order to identify the cell types of the samples in the dataset GSE141445 and describe the cellular states of the clustering results, annotations were made for the seven different cell clustering results. Among the annotated cells, the engagements between endothelial cells and epithelial cells were the most frequent ( Figure 9A , Supplementary Figure 6A ). Epithelial cells and myeloid cells had the strongest interactions with other cells ( Figure 9B ). Epithelial cells had the strongest interaction with endothelial cells, T cells, myeloid cells, and B cells ( Supplementary Figure 6B ). Interaction between epithelial cells and B cells was carried out by MIF−(CD74+CXCR4) ( Figure 9C ). Through reduction and clustering, epithelial cells were eventually categorized into 12 clusters ( Figures 9D, E ). Further investigation into the developmental trajectory of the key cell, epithelial cell, revealed that epithelial cells varied from dark blue to light blue during differentiation and were categorized into 9 stages and 12 clusters ( Figure 9F ). During the differentiation of epithelial cells, only the expression level of FASN continuously increased, and it remained relatively active throughout the entire cell development stage ( Figure 9G ). In conclusion, the development of PCa might be related to epithelial cells and FASN. Furthermore, the expression patterns of prognostic genes during the pseudotime process of myeloid cells differentiation were demonstrated. As shown in the figure, the FASN and TLR3 genes exhibited higher expression levels at the early stages of myeloid cell differentiation, while RACGAP1 and NR3C1 genes had higher expression at the late stages, and BMP2 gene showed higher expression at the mid-stage of differentiation. These results suggested that these five prognostic genes had elevated expression at specific stages of myeloid cell differentiation, which might have led to their overall expression being less prominent ( Supplementary Figure 7 ).

In order to determine the expression patterns of prognostic genes in PCa samples and control samples, expression analysis was conducted using the dataset and RT-qPCR experiments. The expression levels of BMP2, NR3C1 and TLR3 were all remarkably lower in PCa samples than in the control samples (P < 0.0001), while the expression levels of RACGAP1 and FASN were remarkably higher in PCa samples than in the control samples (P < 0.0001) ( Figure 10A ). In RT-qPCR, the expression levels of NR3C1 (P < 0.0001), BMP2 (P < 0.01) and TLR3 (P < 0.01) in PCa were significantly lower than those in control group, while the expression levels of RACGAP1 (P < 0.05) and FASN (P < 0.01) in PCa were significantly higher than those in control group ( Figure 10B ). The expression level and trend of prognostic genes in vitro samples were consistent with that of bioinformatics analysis, indicating the reliability of bioinformatics analysis results. Previous studies have highlighted the significance of BMP2 and FASN in PCa progression. Horvath et al. demonstrated that reduced BMP2 expression is associated with PCa progression, linking its loss to more aggressive phenotypes (52). Similarly, Tae et al. reported that decreased BMP2 expression correlates with a higher incidence of biochemical recurrence (BCR) and elevated Gleason scores (GS) (53). However, the precise mechanisms underlying BMP2’ s role as a tumor outcome determinant remain elusive. In contrast, the biological role of FASN as a key regulator of lipid metabolism is well-established. By driving lipid synthesis, FASN provides energy to fuel tumor proliferation and progression. Despite this, its specific regulatory effects within the PCa microenvironment are not fully understood. Chianese et al. observed FASN overexpression in PCa and proposed that FASN inhibition disrupts the metabolic axis, leading to lipid accumulation and subsequent lipotoxicity (54, 55). This metabolic dysregulation impairs replication mechanisms and arrests cells in the G0/G1 phase, thereby inhibiting proliferation (56, 57). These findings underscore the multifaceted roles of BMP2 and FASN in PCa biology and warrant further investigation into their underlying mechanisms and therapeutic potential.

Western blot analysis revealed that compared with the Vector group, the protein level of BMP2 was significantly upregulated in the OE-BMP2 group, while the FASN protein level showed no obvious change; conversely, in the sh-FASN group, the FASN protein level was notably downregulated compared with the sh-NC group, with no significant difference in BMP2 level ( Figures 11A-D ). Immunofluorescence staining demonstrated that the proportion of Ki67-positive cells was significantly higher in the OE-BMP2 group than in the Vector group ( Figures 11E, G ), and lower in the sh-FASN group than in the sh-NC group ( Figures 11F, H ), as quantified by statistical analysis. In the cell migration assay, the number of migrated cells in the OE-BMP2 group was significantly greater than that in the Vector group at both 24 h and 48 h time points, whereas the sh-FASN group showed a significantly reduced number of migrated cells compared with the sh-NC group at the same time points, as shown by the quantitative results ( Figures 11I-L ).