The PCa data obtained through scRNA-seq was retrieved from the NCBI Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The dataset used for single-cell analysis included 32 tumor and non-tumor samples from 18 PCa patients, with the GSE accession number GSE181294. Detailed information about the samples, including Tumor, Adj-normal, grade, Gleason, Grade Group, Path Stage, Margin, Age, etc., can be found in Supplementary Table S1. Bulk RNA-seq datasets and clinical data were obtained from the Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). The data used in this study was obtained from publicly available databases and therefore did not undergo ethical review.

The scRNA-seq data was imported into R software (version 4.2.0) and examined using the Seurat package (version 4.3.0) (18–21). We performed rigorous quality control on the data excluding the following cells: (1) 300 nFeature < 7,500; (2) 500 nCount < 100,000; (3) mitochondrial gene expression exceeding 25% of the total gene count within the cell; (4) red blood cell gene expression exceeding 5% of the total gene count within the cell. We kept 127,930 cells ultimately for more study. Our research did not call for ethical permission as we made use of publicly accessible database data.

Normalization and selection of the top 2000 highly variable genes (22–26) were performed on the filtered samples using the “NormalizeData” and “FindVariableFeatures” functions in the Seurat package. To correct for batch effects between datasets, principal component analysis was conducted using the Harmony R package (version 0.1.1) (27, 28). Cells were clustered using the FindClusters function with a resolution of 1.0 based on the top 30 principal components (PCs). The top 30 significant PCs were selected for uniform manifold approximation and projection (UMAP) dimensionality reduction and visualization of gene expression (29, 30).

The scRNA-seq data was analyzed for CNVs using the inferCNV R program (version 1.6.0) available from the GitHub repository of the Broad Institute (https://github.com/broadinstitute/inferCNV). This software application facilitates the differentiation between malignant and healthy cells by examining the chromosomal positions and gene expression levels to ascertain CNVs (31, 32). Tumor-EPCs were identified as cells with high CNV scores.

We used Seurat’s “FindAllMarkers” function (33) to conduct a Wilcoxon rank-sum test with the goal of identifying differentially expressed genes (DEGs) among various cell clusters in order to examine the heterogeneity of PCa cells (34). Threshold = 0.25, min.pct = 0.25, and min.diff.pct = 0.25 were the parameters that were employed. Then, we used Seurat’s “DotPlot” and “featureplot” programs to show the expression patterns of the DEGs in each cluster. Manually reviewing the outcomes and consulting pertinent literature helped with the cell annotation process. Additionally, we regrouped these cells in order to investigate the heterogeneity of malignant cells in greater detail. We used marker identification to characterize each subtype based on the genes unique to that grouping.

The Gene Ontology (GO) (35, 36) and Gene Set Enrichment Analysis (GSEA) (37) tools were used to conduct enrichment studies of DEGs in distinct cell types. We performed a functional analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) utilizing the ClusterProfiler R package (version 4.6.2) (38–42).

One technique for gene set enrichment analysis is gene set variation analysis (GSVA). In order to determine enrichment scores for each gene set in each sample, it evaluates the variability of gene expression data and compares it to predefined gene sets.

In order to examine the variations in development and differentiation amongst malignant subtypes in PCa, we created pseudotime trajectories of the subtypes using Monocle (v2.24.0), which showed the patterns of malignant cell differentiation (43). We utilized the CytoTRACE R package (version 0.3.3) (44) to evaluate cell fate, which enabled us to deduce the time course of cell differentiation. Next, we employed the Slingshot software (version 2.6.0) to deduce cell lineages as malignant subtypes differentiated (45). Each malignant subtype’s differentiation trajectory was deduced using the “getlineage” and “getCurves” tools. After that, a UMAP projection was created using these trajectories for visualization.

We used the CellChat R package (version 1.6.1) (46) to compute regulatory networks based on ligand-receptor levels and infer complicated cell-to-cell interactions, as well as analyze the intercellular communication network within the TME. The “Identify Communication Patterns” function was used to estimate the number of communication patterns, and the “netVisualDiffInteraction” function was used to show the variation in communication strength between cells. The P-value, or significance threshold of 0.05, was used.

To analyze the scRNA-seq data and reconstruct gene regulatory networks, we used Python (v3.7) and the pySCENIC package (version 0.10.0) uncovering key gene regulatory mechanisms. To evaluate TFs enrichment and regulatory factor activity, we created an AUCell matrix for this investigation.

To construct a risk model, we obtained PCa-related data from the TCGA database (https://portal.gdc.cancer.gov/). Initially, we performed univariate Cox regression analysis to screen for potential prognostic-related genes (47–52). To account for multicollinearity among these genes, we further employed the least absolute shrinkage and selection operator (LASSO) regression (53–57) using the glmnet package (version 4.1-6). We assigned weights to each gene based on its expression levels and the coefficients obtained from the multivariable Cox regression analysis. This allowed us to construct a risk score formula: Risk score = ∑_i^n (Xi × Yi) (X: coefficient, Y: gene expression level). Next, using the “surv_cutpoint” function, we determined the optimal cutoff value to divide patients into high and low-risk groups based on the risk scores. The predictive accuracy of the model was assessed using ROC curve analysis (58). To observe the prognosis of patients in different groups, Kaplan-Meier survival curves were used to evaluate the survival differences between different risk groups (58, 59). Survival analysis was performed using the Survive package (version 3.3.1) and survminer package (version 0.4.9). To validate the predictive capabilities of the model based on NEFH+ malignant cell scores, we employed the “Survival” and “Time ROC” R packages to generate ROC curves for 1-year, 3-year, and 5-year survival rates, and calculated the area under curve (AUC) value of the ROC curve (60–65). Survival analysis and time-dependent ROC analysis were used for model validation. The distribution of risk score scores, scatter plots of survival status, and heatmaps were utilized to assess the model.

For the evaluation of the TME, we employed the ESTIMATE R package (version 1.0.13) to estimate the stromal score, immune score, and ESTIMATE score in PCa tissues (66). To analyze RNA-Seq data and determine the relative proportion of infiltrating immune cells, we utilized the Cell Type Identification for Estimating Relative Subtypes of RNA Transcripts (CIBERSORT R package, version 0.1.0) algorithm (67, 68), which provides insights into 22 different immune cell types. Additionally, to quantify immune cell infiltration in each sample, we employed the xCell package to assess the enrichment of immune cells in PCa samples. Furthermore, we investigated the correlation between risk scores and immunomodulatory genes, particularly immune checkpoints. To evaluate the response to tumor immune therapy, we utilized the Tumor Immune Dysfunction and Exclusion (TIDE) tool (http://tide.dfci.harvard.edu).

Using the TCGAbiollinks and maftools R package, the somatic mutation data of PCa were obtained from the TCGA database. The PCa expression data and TMB file were imported into the R package. The correlation between DEGs associated with PCa and tumor mutation burden was computed. Based on the correlation results, waterfall plots were generated to depict the high-risk and low-risk groups.

With the pRRophetic R package (version 0.5), we employed the GDSC database (https://www.cancerrxgene.org/), the most comprehensive pharmacogenomics database, to predict the treatment response of each tumor sample (69, 70). The IC50 values for each medication were determined using regression, and the accuracy of both regression and prediction was assessed using 10-fold cross-validation with the GDSC training set. All settings were set to their default values, including the “combat” option for removing batch effects and the mean value for duplication gene expression (71–73).

The MDA PCa 2b and VCap cell lines were acquired from ATCC and cultured in specific growth media. MDA PCa 2b cells were cultured in F12K medium, while VCap cells were cultured in DMEM medium. Both media were supplemented with 10% fetal bovine serum from Gibco BRL, USA, and 1% penicillin/streptomycin. The cells were incubated under standardized conditions, including a constant temperature of 37°C, a carbon dioxide concentration of 5%, and a humidity level of 95%.

IRX4 knockdown was achieved by employing short interfering RNA constructs obtained. The transfection process followed the prescribed guidelines supplied by Lipofectamine 3000 RNAiMAX (Invitrogen, USA). The cells were introduced to a six-well plate when they reached a 50% coverage and subsequently underwent transfection using negative controls (si-NC) as well as knockdown constructs (Si-IRX4-1 and Si-IRX4-2). The transfection process utilized Lipofectamine 3000 RNAiMAX (Invitrogen, USA) for each instance. The siRNA sequences are recorded in Supplementary Table S2 in the Supplementary Materials.

The CCK-8 assay was used to evaluate the cellular vitality of MDA PCa 2b and VCap cells after transfection. The cell suspensions were added to a 96-well plate with a density of 5×10³ cells per well and incubated for 24 hours. After adding 10μL of CCK-8 labeling reagent (A311-01, Vazyme) to each well, the plate was placed in a light-protected environment and incubated at 37°C for 2 hours. Cell viability was assessed by quantifying the absorbance at 450nm using an enzyme-linked immunosorbent assay (ELISA) reader on days 1, 2, 3, and 4. The mean optical density (OD) values were computed and depicted on a line graph.

MDA PCa 2b and VCap cells that had been genetically modified were placed in a 6-well plate at a concentration of 5×10³ cells per well and let to grow overnight. Afterwards, a solution of EdU with a concentration that is twice as strong as the original was prepared by combining 10 mM EdU with a medium that does not contain serum. The solution was introduced into the cell culture and permitted to incubate at a temperature of 37°C for a duration of 2 hours. After the incubation period, the liquid containing the cells was removed, and the cells were carefully rinsed with PBS. Subsequently, the cells were treated with a 4% paraformaldehyde solution for a duration of 30 minutes to ensure fixation. Afterwards, a glycine solution with a concentration of 2 mg/mL and 0.5% Triton X-100 was administered for a duration of 15 minutes. Subsequently, the cells were incubated at room temperature with a solution containing 1 ml of 1X Apollo and 1 ml of 1X Hoechst 33342 for a duration of 30 minutes. Fluorescence microscopy was used to measure and analyze cell proliferation.

The cells were transfected, then grown to a cell density of 95% by seeding them onto 6-well plates. Then, using a sterile 200 μL pipette tip, a careful linear scratch was made over the cell layer in the culture wells. After that, PBS was used to gently wash the wells. After rinsing, the cell culture was let to continue and the culture medium was changed. Photography was used to record the scratches at the original time point (0 hours) and 48 hours later. Measurements were made of the scratches’ width for further examination.

To prepare for the experiment, the cells were subjected to a 24-hour serum-free medium starvation. Subsequently, cell suspensions were added to the upper chamber, which contained Costar, after the addition of matrix gel (BD Biosciences, USA). In the lower chamber, serum-containing medium was added. The cells were then incubated in a cell culture incubator for 48 hours. After the incubation period, the cells were fixed with 4% paraformaldehyde and stained with crystal violet to evaluate their invasive capacity.

We performed statistical analysis using R software and Python software to analyze the database data. All p-values reported in this study are two-tailed, with values less than 0.05 considered statistically significant. P-values below 0.001 were considered highly significant, while those below 0.0001 were regarded as extremely significant.

We conducted an extensive analysis of the obtained dataset to unveil the intricate single-cell landscape within the PCa microenvironment. Our workflow was illustrated in Figure 1. Using the known typical cell type marker genes, 127,930 high-quality cells were labeled. By using dimensionality reduction clustering with UMAP plot, 9 cell types were obtained: T-NK cells, B-plasma cells, endothelial cells (ECs), myeloid-cells, mast cells (MCs), epithelial cells (EPCs), pericytes-smooth muscle cells (pericytes-SMCs), fibroblasts, and plasmacytoid dendritic cells (pDCs). The proportions of different tissue types (Adjacent-Normal high-grade (NHG), Adjacent-Normal low-grade (NLG), Tumor high-grade (THG, Gleason 8–10), Tumor low-grade (TLG, Gleason 6 and 7) and cell cycle phases (G1, G2/M, S) in various cell types were visualized using pie charts (Figures 2A, B).

Subsequently, we analyzed the proportion of cell types in different phases (Figures 2C, D). We found that EPCs, pericytes-SMCs, fibroblasts, ECs, and myeloid-cells were mainly derived from tumor cells. The results of cell cycle study showed that T-NK cells occupied a larger proportion in G2/M and S phases. Then, we described the expression level and distribution of typical marker genes related to cell subtypes in cells (Figures 2E, F). By visualizing the analysis of G2/M.score, S.score, nFeature-RNA and nCount-RNA of all cells, the differences among cell types were further clarified (Figure 2G). The results of enrichment analysis showed that the marker genes of EPCs were mainly enriched in oxidative phosphorylation pathway (Figure 2H). It is worth noting that PCa is mainly transformed from glandular epithelial cells of prostate (74), and oxidative phosphorylation is considered as a tumor-related metabolic marker (75), which confirms that our findings are consistent with the recognized biological functions related to PCa.

Given the profound importance of malignant cells in TME, our subsequent objective is to characterize these cells in the microenvironment of PCa. Studies indicated that a high CNV score typically signified significant copy number alterations within cells, which might have been associated with malignant transformation. CNV scoring could be utilized to identify tumor heterogeneity and reveal the presence of different subpopulations within the tumor microenvironment, aiding in the understanding of tumor development and resistance mechanisms (76, 77). To detect aberrant amplification or deletion of chromosome copy number in EPCs, we initially employed inferCNV to analyze the chromosome CNV of epithelial cells using endothelial cells as a reference (Supplementary Figure S1). According to CNV level, malignant cells was distinguished from EPCs. After that, we re-clustered 3,076 malignant cells, and annotated them according to each cell marker gene, and identified four malignant cell subtypes: C0 TRPM4+ malignant cells, C1 SLPI+ malignant cells, C2 MIR205HG+ malignant cells and C3 NEFH+ malignant cells (Figure 3A). We employed the UMAP plot combined with a cellular proportion pie chart to illustrate the relative distribution of different subgroups (THG and TLG) within the four malignant cell subtypes (Figure 3B). Additionally, we visualized these cells based on tissue types (Figure 3C). The results demonstrated that, compared to other subtypes, the C1 and C3 malignant cell subtypes exhibited a higher proportion of THG tissue. Additionally, a stacked bar graph revealed that THG tissue predominantly comprised the C1 and C3 subtypes, in contrast to TLG tissue (Figures 3D, E). Therefore, we hypothesized that the heterogeneity between these two tissue types may be associated with the C1 and C3 subtypes. Similarly, the Ro/e preference plot indicated a higher cell abundance of C1 and C3 subtypes in THG tissue, further substantiating our conclusion (Figure 3F). Figure 3G showcased the differential expression of the top five marker genes in the malignant cell subtypes, visualized using the volcano plots (Figure 3H).

Next, we showed the results of CNVscore, nFeature-RNA and nCount-RNA of different tissue types and malignant cell subtypes by violin plots (Figures 3I, J). It was found that the CNVscore of THG was higher than that of TLH, which was consistent with the biological process of PCa tissue. It is noteworthy that both the C1 and C3 subtypes exhibited higher CNV scores. Additionally, compared to the C1 subtype, the C3 subtype with high NEFH expression showed elevated levels of nCount-RNA and nFeature-RNA. Therefore, we inferred that the C3 subtype’s malignancy level might have been higher.

Subsequently, we determined that the subtypes shared specific biological functions through the enrichment analysis of a variety of malignant cells. For instance, transition metal homeostasis was associated with both C1 and C3 subtypes. Furthermore, the C3 subtype was enriched in the cholesterol metabolism process (Figures 3K–M). We also illustrated the metabolism pathways using bubble plot and discovered that the C3 subtype was substantially enriched in the metabolic pathways of glutathione metabolism and fatty acid biosynthesis (Figure 3N).

To understand the source and development of cancer cells, we extensively studied the intricate lineage and advancement of malignant cells utilizing the Monocle software. It was easy to see in Figures 4A–E that the C2 MIR205HG+ malignant cells were mostly in the early stages of differentiation, more specifically in state 1 along the time axis. Conversely, the C3 NEFH+ malignant cells were in the last stage of development, predominantly found within state 3 throughout the temporal trajectory. Afterwards, we utilized the CytoTRACE technique to evaluate the differentiation and developmental correlation between several subtypes of tumor cells. The findings indicated that C3 NEFH+ malignant cells displayed elevated cellular stemness, as seen in Figure 4F. Combining the Slingshot analysis, we discovered a differentiation boundary transitioning from TLG to THG tissue types, indicating a higher malignant degree among the cells at the terminal stage of differentiation (Figure 4G). Malignant cells often possess self-renewal capability and differentiation potential. Thus, as the tumor advances, malignant cells in the last stage of differentiation tend to have greater cellular stemness, in line with the results obtained from the CytoTRACE investigation. Furthermore, the progression of cancerous cell subtypes can be described in the order of C2→C1→C0→C3. Additionally, the slingshot analysis of time states showed that State 3 was positioned at the end of one of the differentiation branches, with a greater proportion of C3 subtypes within state 3. This finding further supports the conclusions drawn from the Monocle analysis (Figures 4H, I).

Consistent with those findings, the genes TRPM4, which serve as markers for the C0 subtype, and NEFH, which serve as markers for the C3 subtype, were predominantly expressed during the mid-late stage of the developmental trajectory. On the other hand, the expression levels of SLPI, a marker for the C1 subtype, and MIR205HG+, a marker for the C2 subtype, were initially high but declined over time (Figure 4J). In addition, we conducted GO-BP enrichment analysis on the DEGs linked to the malignant cell subtypes. The analysis revealed that these genes were mainly enriched in biological processes related to immunity, cytotoxicity, antigen, processing, and other similar activities (Figure 4K).

We utilized CellChat to infer and analyze communication between tumor cell subtypes and other cell types from single-cell data (Supplementary Table S3). The number and intensity of interactions between all cell types in PCa samples were comprehensively summarized (Figure 5A). It was found that compared with other types of cells, C3 NEFH+ malignant cells had a more significant effect on pericytes-SMCs and fibroblasts. The circle graphs quantified the number and intensity of interactions between all cells with C3 NEFH+ malignant cells as the signal source and fibroblasts as the target respectively (Figures 5B, C). The results showed that there was a strong intercellular communication network between C3 NEFH+ malignant cells and fibroblasts.

Next, we identified the ligand-receptor signals associated with the communication pathway (Figure 5D) to determine the primary afferent and efferent signals related to the C3 NEFH+ malignant cell subtype and other cells. The findings indicated that the primary ligands associated with the output of C3 NEFH+ malignant cells, when employed as signal senders, were MIF, MK, GDF, and CD46. As signal receivers, the fibroblast-related receptors mainly included PTN, CD99, and PDGF.

Subsequent analysis revealed potential connections to the PTN signaling pathway network. Through network centrality analysis of the inferred PTN signaling network, we found that C3 NEFH+ malignant cells can act as signal senders within the PTN pathway. Fibroblasts, on the other hand, can function both as signal sender promoting their transformation into cancer-associated fibroblasts (CAFs) and as signal receiver, mediator, and influencer interacting with C3 NEFH+ malignant cells (Figure 5E). Notably, C3 NEFH+ malignant cells demonstrated the ability to engage in paracrine interactions with fibroblasts, resulting in a substantial communication intensity between these cell populations (Figure 5F). In addition, we compared the receptor-ligand interaction between C3 NEFH+ malignant cells and other cell types and found that when this subtype interacted with fibroblasts, the ligand receptor had a high communication probability with PTN-NCL (Figures 5G, H). Additionally, a circle graph further confirmed that the interactions between C3 NEFH+ malignant cells and fibroblasts could be mediated through the receptor-ligand pairs within the PTN signaling pathway, specifically involving PTN-NCL (Figure 5I).

Essentially, our study provided profound insights into the intricate interactions between fibroblasts and malignant cell subtypes in PCa. This relationship is likely closely linked to the transformation of fibroblasts into CAFs, which promotes the progression of PCa.

TFs can directly interact with the genome and regulate gene transcription by binding to specific nucleotide sequences upstream of the target gene. This interaction plays a significant role in determining the biological functions of cells.

To begin with, we employed the SCENIC and connection specificity index matrix to classify prostate malignant cells into four regulatory modules (M1, M2, M3, M4) based on the similarity of AUCell score rules (Figure 6A). Subsequently, we conducted dimensionality reduction and clustering analyses considering various subtypes and tissue types (Figures 6B, C).

Through a comparison of the expression levels and regulatory activities of TFs within each module and the malignant cell subtypes, we identified that the TFs in the M4 module predominantly regulated C3 NEFH+ malignant cells (Figures 6D, E). Within the M4 module, we found that IRX4 exhibited the highest fraction of variance across the subtypes, indicating its prominent role in explaining a significant portion of the data variability within the M4 module (Figure 6F). Importantly, IRX4 demonstrated a high specificity score in C3 NEFH+ malignant cells and THG (Figures 6G, H). This suggested a strong and specific regulatory relationship between IRX4 and its target genes, highlighting its potential as a biomarker or therapeutic target.

Finally, we visualized the expression levels of five key regulatory factors (NEUROD1, IRX4, LTF, PURA, and ELK4) in the C3 subtype (Figures 6I, J). We observed that the expression of IRX4 in the C3 subtype was significantly higher compared to other malignant cell subtypes. To further investigate the relevance of these TFs, we performed survival analysis using Kaplan-Meier and AUC curves (Supplementary Figures S2A–F). Interestingly, our results indicated that high expression of IRX4 might be associated with a poorer prognosis in PCa. Nevertheless, the specific mechanism by which IRX4 influences PCa remains unclear. Therefore, conducting in vitro functional experiments to validate the impact of IRX4 on PCa cells is imperative.

To further investigate the role of IRX4 in PCa, we conducted in vitro experiments using MDA PCA 2b and VCap cell lines. Initially, we knocked down IRX4 and assessed the mRNA and protein expression levels before and after knockdown. Our findings revealed a significant decrease in mRNA and protein expression levels in both cell lines compared to the control group (Figure 7A). Moreover, there was a noticeable decrease in cell viability following the knockdown (Figure 7B).

Subsequent colony experiments demonstrated a significant reduction in the number of cells after IRX4 knockdown (Figures 7C, D). Additional EDU experiment confirmed that the knockout of IRX4 partially inhibited cell proliferation (Figure 7E). Furthermore, the wound healing assay and Transwell assay indicated a substantial decrease in cell migration after IRX4 knockdown (Figures 7F–H), and the cell invasion capability also decreased (Figure 7I).

Collectively, these results indicate that the knockdown of IRX4 can inhibit the activity, migration, and proliferation of tumor cells, thereby impeding tumor growth.

We developed a prognostic model to investigate the clinical significance of the NEFH+/IRX4 regulatory network. Initially, we performed univariate Cox regression analysis to identify genes significantly associated with prognosis (Figure 8A). To address the issue of multicollinearity among these genes, we employed LASSO regression analysis for further selection (Figure 8B). Subsequently, a multivariate Cox regression analysis was conducted, resulting in the identification of five genes related to prognosis. The coef values for these genes were calculated (Figures 8C, D). The findings revealed that ZNF782, ZNF695, YY1, NR1I3, and FOXA3 were unfavorable prognostic factors.

To further investigate the differences between different scoring groups, we performed an analysis of DEGs. Based on the optimal cut-off value of NEFH+ malignant cell score, patients in the TCGA cohort were categorized into two groups: the high NmRS group and the low NmRS group (NmRS: NEFH+ malignant cells risk score). It was observed that higher scores were associated with worse prognosis. Curve and scatter plots were utilized to illustrate the differences in risk scores, survival, and outcomes between the two groups, clearly indicating that the high NmRS group was associated with a poorer prognosis (Figure 8E). Furthermore, a heatmap was generated to display the differential expression of the five genes between the high and low NmRS groups (Figure 8F). Principal component analysis demonstrated that PC1 (high NmRS group) accounted for 10.2% of the total variance in all principal components, while PC2 (low NmRS group) accounted for 4.2% of the total variance (Figure 8G). Additionally, the ROC curve provided an intuitive visualization of the AUC values predicted by the TCGA cohort at 1 year, 3 years, and 5 years, demonstrating the predictive value of the model (Figure 8H). The Kaplan-Meier survival curve further confirmed the conclusion that the high NmRS group had a worse survival outcome, with a p-value less than 0.0001 (Figure 8I).

To elucidate differential gene expression and associated biological processes between high and low groups, we employed visualization and enrichment analysis techniques. Initially, a heatmap was utilized to depict gene expression of the top thirty genes (Figure 9A), while a volcano plot showcased differential gene up-regulation and down-regulation (Figure 9B).

Subsequently, various enrichment methods were employed to delve deeper into the related biological processes. KEGG analysis indicated predominant enrichment of differential genes in cytoskeleton in muscle cells, motor proteins, and cardiac muscle contraction (Figure 9C). Conversely, GO analysis revealed enrichment in myofibril and muscle cell development (Figure 9D). Furthermore, GSVA enrichment analysis was conducted on the gene set comprising the prediction model, as represented in the heatmap (Figure 9E). Lastly, GSEA analysis was performed on the initial 30 up-regulated and down-regulated genes. Up-regulated genes exhibited enrichment in DNA geometric change, sister chromatid segregation, mitotic sister chromatid segregation, cell cycle DNA replication, down-regulated genes were primarily enriched in striated muscle adaptation, sarcomere organization, striated muscle cell development, and muscle cell development (Figure 9F).

To investigate immune cell composition differences in NmRS with varying risk scores, we employed the CIBERSORT algorithm for analyzing twenty-two immune cell types in TCGA database PCa patients, as illustrated in Figure 10A. The proportions of eight immune cell categories were presented in Figure 10B, revealing distinctions between the two groups via a box plot. Notably, higher infiltration of naive B cells, T cells CD4 memory resting, and T cells follicular helper was observed in the high NmRS group, implying potential immune reactions associated with immune escape.

Subsequently, we assessed the correlation between immune cells and NmRS, as depicted in Figure 10C. The findings indicated a significant positive correlation between NmRS and macrophages, T cells CD4 memory resting, and T cells follicular helper, while a significant negative correlation was observed between NmRS and M1 macrophages. A heatmap visualization was employed to depict the correlation analysis among immune cells, modeling genes, OS, and risk score (Figure 10D). Moreover, differences in TIDE values between the two groups were evident (Figure 10E).

Additionally, Figure 10F displayed a heatmap showcasing variations in modeling genes, StromalScore, ImmuneScore, ESTIMATEScore, TumorPurity, and immune cell infiltration levels between the high NmRS and low NmRS groups, computed using the CIBERSORT and Xcell algorithms.

Next, we presented the correlation between DEGs and tumor mutation load using a waterfall diagram for the high and low-risk groups (Figure 10G). Furthermore, the correlation between immune checkpoint-related genes, modeling genes, risk score, and OS was displayed in a bubble plot (Figure 10H). The results revealed a strong positive correlation between NR1I3, ZNF782, YY1, and most immune checkpoints.

Finally, through drug sensitivity analysis, we identified potential clinical efficacy of certain drugs for prognosis-related genes, including ABT.263, AZD.2281, Bleomycin, Cisplatin, Cytarabine, Embelin, Epothilone.B, Etoposide, IPA.3, Methotrexate, MS.275, Shikonin, and Vorinostat. However, for the low-risk group, Dasatinib and Bicalutamide were found to be effective (Figure 10I).