The 13 PCa samples with scRNA-seq data (GSE141445) were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih/) and were used to identify MRMGs. In addition, a total of 548 samples were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) including 496 PCa samples and 52 normal samples. After eliminating incomplete data, 346 PCa samples were included in the analysis as a training cohort. Subsequently, to evaluate the predictive accuracy of the model in a verification cohort, GSE70768 (n=110) and GSE46602 (n=36) collected from the GEO database were merged and eliminated batch effects by R package “inSilicoMerging” and empirical Bayes methods (23).

We used the “Seurat” package in R to generate objects and filter out high-quality cells. The filtering criteria were to remove the genes from less than 3 cells, the cells with less than 50 genes detected, or the cells with more than 5% mitochondria genes. The data were then normalized using the “LogNormalize” method, and the first 1500 highly variable genes identified by JackStraw analysis were analyzed by principal component analysis (PCA). Furthermore, the “FindClusters” function (resolution=0.5) was used to cluster and visualize the first 15 PCA in t-distributed stochastic neighbor embedding (t-SNE). Finally, the marker genes of each cluster were identified (cut-off thresholds: adjusted P-value < 0.05 and |logFC| > 1) by comparing the differences of gene expression between a cluster and all other clusters using the “FindAllMarkers” function and Wilcoxon-Mann-Whitney test, and the cell types were annotated and visualized by the “SingleR” package.

Based on the MRMGs, we utilized the “ConensusClusterPlus” package in R for consensus molecular clustering methods and divided the patients from the TCGA cohort into different clusters. After 1000 initial resamples and 50 iterations by the unsupervised clustering algorithm, the most optimal clustering number was selected according to the consistency matrix, the cumulative distribution function (CDF), and the relative change of the area under the CDF curve.

Subsequently, the differences and bRFS rates among different clusters were computed and visualized by the R packages of “scatterplot3d”, “survival”, and “survminer”. The relationships between different clusters and the clinicopathological characteristics (age, Gleason score, prostate specific antigen (PSA), pathological T stage, and pathological N stage) were analyzed with the chi-square test. Besides, the differences of somatic mutations from TCGA in different clusters were analyzed using the R package “maftools”, and a single-sample GSEA (ssGSEA) algorithm with the “GSVA” package in R was utilized to estimate the infiltration of immune cells in each sample based on the gene sets of immune cell types obtained from the previous study (24).

Based on the differentially expressed genes (DEGs) identified via the “LIMMA” package in R (FDR < 0.05) among three MRMGs clusters or between high- and low-risk groups, the “clusterProfiler” R package was used to perform gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (25), and P-value < 0.05 was considered to be significant enrichment. P-value was adjusted with the Benjamini-Hochberg methods.

Gene set enrichment analysis (GSEA) was used to analyze the relative richness of specific gene sets in the sample population, and to detect statistical differential expression patterns among three MRMGs clusters or between high- and low-risk groups, with FDR <0.05.

Based on the MRMGs identified by scRNA-seq analysis, the prognostic MRMGs were identified as candidate genes by univariate Cox regression analysis. The DEGs between PCa tumors and normal samples from the TCGA were selected through Wilcoxon test by the R package of “LIMMA” (FDR < 0.05). Then the intersection of scRNA-MRMGs, prognostic MRMGs, and DEGs was obtained using the “Venn” package in R. Meanwhile, a protein-protein interaction (PPI) network of the intersection was built by the STRING database (https://string-db.org/) (26) and visualized by Cytoscape (version 3.8.2) (27). The intersection genes were evaluated by least absolute shrinkage and selection operator (LASSO) Cox proportional hazards regression using the “glmnet” package to eliminate the effect of overfitting. Then, the criterion for identifying optimal model genes: According to the minimum Akaike information criterion (AIC), the optimal model genes were obtained by taking the optimal penalty parameter (λ) corresponding to the minimum 1-standard error (SE) obtained from the 10-fold cross-validation results (28).

Subsequently, the risk score was calculated by summing the expression and coefficient of each gene of MRMGPS. The distribution of signature genes within different cell types was further visualized by the R packages of “Vlnplot”, “Dimplot” and “Featureplot” at the single-cell level, and the CIBERSORT algorithm was used to estimate the relative proportion of macrophage proportions deconvoluted from the TCGA cohort and then calculate the correlation with gene expression. According to the median value of risk score, patients were divided into two groups with low- or high-risk score. Furthermore, t-SNE and PCA were performed to reveal distribution between two groups by the packages of “stats” and “Rtsne”. Receiver operating characteristic (ROC) curves and Kaplan-Meier curves were generated to evaluate the accuracy of MRMGPS by “timeROC”, “surviminer” and “survival” in R. To further verify the accuracy and reliability of the novel signature, the merged dataset (GSE70768 and GSE46602), named GEO-Merged cohort, was conducted to verify the reliability of the signature in the same way as the TCGA cohort.

20 fresh PCa and 20 benign prostatic hyperplasia (BPH) tissue samples were collected from the Third Affiliated Hospital of Sun Yat-sen University (Guangzhou, China) from 2021 to 2022. All the selected samples were informed consent of the patients and approved by the Clinical Ethics Board of the Third Affiliated Hospital of Sun Yat-sen University. Following a standard protocol, total RNA was extracted by TRIzol reagent (Invitrogen, Cat #15596018) and reverse transcribed by SuperScript™ III Reverse Transcriptase (Invitrogen, Cat # 18080044). SYBR Green qPCR reagent (GenStar, Cat #A311) was applied for RT-qPCR. GAPDH was used as the internal control. The primers used for RT-qPCR were listed in Table S1 .

The IHC data and images collected from the Human Protein Atlas (HPA, https://www.proteinatlas.org/) (29) were used to reveal expression levels of protein of signature genes between PCa and normal samples.

The differences between MRMGPS-based risk score and clinicopathological features were analyzed by a chi-square test, and stratified survival analysis of patients in two subgroups was conducted to evaluate the robustness of MRMGPS in both TCGA and GEO-Merged cohorts.

Meanwhile, independent prognostic indicators were determined by COX proportional hazard regression model and the independent indicators were used to construct nomogram for predicting 1-, 3-, and 5-year bRFS by using the “rms” package. Furthermore, decision curve analysis (DCA), ROC curves, and calibration curves were performed to assess the accuracy and discrimination of the nomogram.

Firstly, the immune infiltration and function of two groups based on MRMGPS were calculated by the EPIC (30), XCell (31), MCPCOUNTER (32), QUANTISEQ (33), CIBERSORT-ABS, CIBERSORT (34), and TIMER (35) algorithms. Secondly, the ssGSEA algorithm with the “GSVA” package and Mann-Whitney test were utilized to calculate the infiltration scores of 29 immune cells or pathways obtained from the previous study (36). The gene sets with annotation were shown in Table S2 . Thirdly, the ESTIMATE scores, immune scores, stromal scores, and tumor purity of the two groups were estimated and compared by the R package of “estimate” (37). Finally, the scores of 7 steps in the cancer-immunity cycle based on MRMGPS were analyzed by Tracking Tumor Immunophenotype (TIP, https://biocc.hrbmu.edu.cn/TIP/) as in previous studies (38).

The Copy Number Variation (CNV) data of all PCa samples in TCGA processed by GISTIC software (39) were downloaded from GDC (https://portal.gdc.cancer.gov/), and further integrated for analysis and visualization. The somatic mutation data in TCGA processed by the VarScan platform were analyzed between two groups by the package of “maftools” and the tumor mutational burden (TMB) between the two groups was compared. In addition, the markers of genomic heterogeneities, including homologous recombination defect (HRD), microsatellite instability (MSI), mutant-allele tumor heterogeneity (MATH), loss of heterozygosity (LOH), ploidy, neoantigens, and RNA-modified genes ( Table S3 ), were collected from previous studies (40). The markers of tumor stemness, including DNA-methylation-based stemness scores (DNAss) and RNA-based stemness scores (RNAss), were collected from the study of Malta TM (41). Furthermore, the correlation between them and MRMGPS-based risk score were analyzed by Pearson or Mann–Whitney methods.

The differences in risk score and 46 immune checkpoints obtained from the previous research (42) were examined to explore the relationship. In addition, the immunophenoscore (IPS) of PCa samples in TCGA were collected from The Cancer Immunome Atlas (https://tcia.at/) (24) and were used to compare the responses of PD1 and CTLA4 in two subgroups.

According to the IC50 value from the Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/) (43), the efficacy of common anticancer drugs between two subgroups was calculated by the “pRRophetic” in R. Besides, the pharmacological information about the anticancer drug targets was collected from DrugBank database (https://go.drugbank.com/) (44), and was used to reveal their expression level in different groups.

The molecular structures of the targeted drugs were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/) (45) and the coordinates of MRMGPS-based genes were obtained from Protein Data Bank (PDB, http://www.rcsb.org/pdb/home/home.do) (46). A Molecular Operating Environment (MOE, version 2019.0102) was used to exclude all water molecules from the protein and ligand files and then add polar hydrogen atoms. Finally, molecular computing and docking were performed to assess the binding energy and interaction mode between targeted drugs and MRMGPS-based genes.

This study used R software (version 4.2.1) for data analysis and graphic visualization. For quantitative data, an independent t-test was utilized to analyze normal distribution variables, and the Wilcoxon test was used to analyze non-normal distribution variables. For qualitative data, a one-way analysis of variance was used to analyze ordinal categories variables, and the Kruskal-Wallis test was utilized to analyze unordered categorical variables. Kaplan-Meier survival analysis was performed by the log-rank test. For all tests, the value of P less than 0.05 was considered statistically significant (*P <0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. ns, not significant).

By analyzing the scRNA-seq data from GSE141445, we collected gene expression profiles of 33441 single cells from 13 PCa samples. Then 21 cell clusters were identified by PCA analysis, and the cells in cluster 10 were defined as macrophage ( Table S5 ). Moreover, a total of 307 genes were differentially expressed in the macrophage cluster, which were defined as MRMGs and used for further research ( Figures 2A–C ). The correlation of 307 MRMGs was presented in the heatmap ( Figure 2D ).

To explore the expression features of MRMGs, we first performed PCA analysis based on the expression of MRMGs between PCa and normal samples in TCGA. The results suggested that individuals could be distinguished from PCa to normal ( Supplementary Figure 1A ), indicating different regulatory roles of macrophages in normal and PCa. Based on the expression of MRMGs, we clustered the consensus of the patients in the TCGA cohort, and according to the area under the line of the CDF curve, the downward trend of CDF delta, and average consistency within the cluster groups ( Supplementary Figures 1B–D ), we obtained optimal three clusters (k=3, C1 = 121, C2 = 99, C3 = 126) ( Figure 3A ). Further PCA analysis showed a distinct clustering, indicating a notable differential distribution of the MRMGs in three clusters ( Figure 3B ). Survival analysis revealed that the patients in C3 suffered the worst prognosis, whereas the patients in C1 had the best prognosis (P < 0.001) ( Figure 3C ). In addition, the analysis of the relationship among different clusters and clinical features showed that C3 had a higher Gleason score (P < 0.001) and an advanced pathological N stage (P = 0.031) ( Figure 3D ), which explained the poor prognosis.

Meanwhile, different clusters also showed remarkable diversity in somatic mutations and immune phenotypes ( Supplementary Figure 1E ). The level of immune infiltration in 13 out of 23 immune cells in C3 was significantly higher than that in C1 and C2 ( Figure 3E ). Hence, we hypothesized that macrophages played a significant role in the features of immunological infiltration based on the aforementioned data.

Subsequently, to explore the differences in biological function among different clusters, we carried out an enrichment analysis based on 4200 DEGs among three MRMGs clusters. These DEGs were not only enriched in biochemical processes related to the regulation of lymphocyte activation and cytokine receptor interaction pathways but also GSEA showed that C3 was likely to be involved in cancer-related pathways explaining the worse prognosis of C3 ( Figures 3F, G ). In general, our data suggested that macrophages played a significant role in the progression of PCa, and we had enough reason to predict that MRMGs might provide important prognostic information.

To select the genes with prognostic value from scRNA-MRMGs for further studying, we first selected 21 overlapping genes among 307 scRNA-MRMGs, 18304 DEGs between PCa and normal samples, and 38 prognostic MRMGs ( Supplementary Figures 2A–D ). Then, through LASSO-Cox regression analysis based on the minimum AIC (λ = 0.03), we finally constructed an optimal signature containing nine genes (CTSZ, FCGRT, GOLM1, SMIM22, ACPP, FAM3B, TFF3, PCA3, and MSMB), named MRMGPS ( Supplementary Figures 2E, F ). Based on their coefficients, MRMGPS-based risk score was built: Risk score = (0.411 × CTSZ) + (0.025 × FCGRT) + (-0.161 × GOLM1) + (-0.107 × SMIM22) + (-0.012 × ACPP) + (-0.029 × FAM3B) + (-0.056 × TFF3) + (-0.025 × PCA3) + (-0.080 × MSMB).

According to the formula for calculating the risk score, it could be considered that the risk score is largely determined by the expression of CTSZ, which was mainly expressed in immune cells such as macrophages and monocytes (47). Therefore, we analyzed the expression profiles of model genes across all cell types identified in the scRNA-seq data and further calculated the correlation between CTSZ expression and macrophage proportions deconvoluted from the bulk RNA-seq data of TCGA cohort. The results showed that CTSZ was mainly expressed in macrophages ( Supplementary Figures 3A–C ) and positively correlated with macrophage proportions (R = 0.11, P = 0.04) ( Supplementary Figure 3D ), suggesting the reliability of our formula for calculating the risk score.

We further examined the mRNA expression profiles of the model genes by RT-qPCR. The result demonstrated that CTSZ, GOLM1, SMIM22, FAM3B, TFF3, and PCA3 had higher expression, while FCGRT, ACPP, and MSMB had lower expression in PCa than BPH samples ( Figure 4A ). We also observed consistent trends in the protein levels of CTSZ, FCGRT, GOLM1, ACPP, FAM3B, TFF3, and MSMB by exploring the HPA database ( Supplementary Figure 4 ). Unfortunately, there were no matching IHC images of PCA3 and SMIM22 in the database. In addition, the network diagram provided a comprehensive view of survival and interactions between model genes ( Figure 4B ).

According to the median risk score, patients were divided into high- or low-risk groups. The risk score plot and survival status indicated that the low-risk group had lower BCR rates ( Figure 4C ). The results of the PCA and t-SNE analysis implied that the distribution of samples in two groups showed two trends ( Figures 4D, E ). Moreover, the Kaplan-Meier curve revealed that patients in the low-risk group had better bRFS (P < 0.001) ( Figure 4F ). Consistently, the area under the ROC curve (AUC) for 1-, 3-, and 5-year bRFS were 0.674, 0.787, and 0.803, respectively ( Figure 4G ). The above results showed that MRMGPS had good accuracy in predicting the prognosis of PCa patients.

To better verify the robustness of MRMGPS, we combined two independent cohorts (GSE70768 and GSE46602) into an external cohort (GEO-Merged) for evaluation ( Supplementary Figure 5A ). After removing the batch effect, the data distribution tended to be consistent between each dataset. The median was homogenized, the mean and variance were comparable, and the samples between the datasets were clustered and intertwined with each other ( Supplementary Figure 5B ), suggesting that the batch effect was better removed.

The same formula as in the TCGA cohort was used to calculate the risk score in the GEO-Merged cohort. The risk score plot and survival status demonstrated consistent trends as shown in the TCGA cohort ( Supplementary Figure 5C ), and indicated that the risk score was a reliable predictor of BCR in patients with PCa. Consistent with the TCGA cohort, PCA and t-SNE also demonstrated excellent separation between two groups ( Supplementary Figures 5D, E ). Besides, survival analysis demonstrated that individuals in the high-risk group had worse bRFS than patients in the low-risk group (P < 0.05) ( Supplementary Figure 5F ). ROC analysis revealed that the AUCs of the MRMGPS were 0.652 at 1 year, 0.685 at 3 years, and 0.700 at 5 years ( Supplementary Figure 5G ).

To explore the association with clinical and pathological characteristics, we first compared the differences in MRMGPS-based risk score according to different stratified characteristics. The violin charts indicated that patients with advanced age (P < 0.05), Gleason score (P < 0.001), pathological T stage (P < 0.01), and pathological N stage (P < 0.001) had higher risk score in TCGA cohort ( Figures 5A–C ). The GEO-Merged cohort showed the same significant differences in Gleason score (P < 0.05) and pathological T stage (P < 0.05) as the TCGA cohort ( Supplementary Figures 6A–C ). Moreover, we found that both TCGA cohort and GEO-Merged cohort patients with advanced Gleason score (TCGA cohort: P < 0.001; GEO-Merged cohort: P = 0.003) and pathological T stage (TCGA cohort: P = 0.003; GEO-Merged cohort: P = 0.008) had poorer bRFS by stratified survival analysis ( Figure 5D , Supplementary Figure 6D ).

Next, to explore the independence of MRMGPS-based risk score and clinicopathological factors, we also carried out univariate and multivariate Cox regression analysis. Univariate Cox regression analysis showed that Gleason score, PSA, pathological T stage, and risk score were significantly correlated with bRFS in both TCGA and GEO-Merged cohorts. Multivariate Cox regression analysis showed that they remained independent predictors of bRFS ( Table S6 ).

To explore the value of MRMGPS in clinical practice, we built a predictive nomogram for bRFS of the independent factors (risk score, pathological T stage, PSA, and Gleason score) obtained by multivariate regression analysis in TCGA and verified by GEO-Merged cohort ( Figure 6A , Supplementary Figure 7A ). As expected, the calibration charts of the nomogram showed excellent predictive performance in both the training and validation cohorts ( Figure 6B , Supplementary Figure 7B ). The AUCs in both cohorts demonstrated that the nomogram was more effective at predicting bRFS than the risk score, PSA, Gleason score, and pathological T stage alone (TCGA cohort: 1-year = 0.709, 3-year = 0.797, 5-year = 0.830; GEO-Merged cohort: 1-year = 0.656, 3-year = 0.759, 5-year = 0.789) ( Figure 6C , Supplementary Figure 7C ). Furthermore, the nomogram had a higher net benefit than other independent factors in two cohorts ( Figure 6D , Supplementary Figure 7D ). The findings indicated that the nomogram for bRFS had excellent predictive ability for long-term BCR risk and brought greater application value.

To investigate potential biological functions and pathways related to MRMGPS in PCa, we screened out 4701 DEGs between two risk groups in the TCGA cohort. GO and KEGG analysis demonstrate that the DEGs were not only enriched in a variety of molecular functions and biological processes related to immunity, such as immune receptor activity, lymphocyte activation, and immune system regulation but also found that the primary immunodeficiency pathway was enriched ( Figures 7A–D ). In addition, GSEA analysis revealed that the high-risk group was mainly enriched in pathways about immune, including primary immunodeficiency, T cell receptor signaling pathway, and NK cell-mediated cytotoxicity ( Figure 7E ). These findings suggested that immunity might be a vital factor in bRFS disparities between subgroups.

To further study the relationship between MRMGPS and tumor immune microenvironment of the PCa samples, we first showed the heatmap and the correlation coefficient of the immune cells with risk score calculated by seven algorithms ( Figures 8A, B ). Additionally, patients in the high-risk group had significantly higher stromal scores, ESTIMATE scores, and immune scores, but the tumor purity showed an inferior to that of the low-risk group (P < 0.001) ( Figure 8C ), indicating that PCa tissue from the high-risk subgroup contained more immune cells and immune molecules. It’s interesting to note that the high-risk group had much higher scores for 13 out of 16 of immune cells and 12 out of 13 of immune-related processes ( Figures 8D, E ). Furthermore, the immune activity scores at each step calculated by TIP analysis revealed that two groups had certain differences in 5 out of 7 steps of the cancer immune cycle ( Figure 8F ). These results revealed the unique tumor immune microenvironment in PCa.

Given tumor heterogeneity has a profound impact on the immune microenvironment, we comprehensively explored crucial features of tumor heterogeneity based on MRMGPS including CNV, mutation profiles, TMB, epigenetic modification, HDR, MSI, MATH, LOH, ploidy, neoantigen, and tumor stemness (40, 41, 48).We found that among MRMGPS-based genes, ACPP, GOLM1, PCA3, SMIM22, and CTSZ showed comparatively high amplification, whereas TFF3, FAM3B, MSMB, and FCGRT showed primarily deletion ( Figures 9A, B ). However, no significant difference in TMB was found between the risk groups (P = 0.24). ( Figure 9C ). The mutation profiles showed that there were more mutations in the high-risk group, and the missense mutations in both groups were higher, in which the mutation rates of TTN and TP53 were both higher than 10% ( Figure 9D ). Furthermore, the correlation between MRMGPS and epigenetic modification-related genes (m1A = 10, m5C = 13, and m6A = 20) was calculated and the result showed that 20 out of 43 genes were significantly related to risk score ( Figure 9E ). Moreover, we observed that risk score revealed a positive correlation with HRD (R = 0.32, P < 0.001), MSI (R = 0.20, P < 0.001), ploidy (R = 0.29, P < 0.001), neoantigen (R = 0.19, P < 0.01) and RNAss (R = 0.23, P < 0.001) but not with MATH and DNAss ( Figures 9F–M ). These results reflected the potential application value of MRMGPS in the exploration of heterogeneity of PCa.

Taking into account the critical role immune checkpoints played in immune responses and therapy, we further examined the relationship between the expression of immune checkpoints-associated genes and risk score. We discovered that 32 of 46 immune checkpoints, including PDCD1 and CTLA4, showed a significant correlation with MRMGPS and 2 checkpoints, HHLA2 and CD44, showed significant negative correlation (P < 0.05) ( Figure 10A ). Additionally, we specifically calculated the IPS of patients undergoing different treatments and showed that patients with low-risk score might be more responsive to CTLA4-positive/PD1-positive immunotherapy (P = 0.032) and had a more favorable therapeutic outcome ( Figure 10B ).

To explore the possibility of applying MRMGPS to the personalized and accurate treatment of PCa, we first compared the differences in sensitivity to different anticancer drugs between high- and low-risk groups according to the IC50 value. The results showed that 10 commonly used anti-PCa drugs showed satisfactory differences in sensitivity, that is, patients in the high-risk group were more likely to be sensitive to bicalutamide and docetaxel, whereas individuals in the low-risk group were more likely to be sensitive to doxorubicin, etoposide, gemcitabine, methotrexate, mitomycin C, paclitaxel, vinblastine, and cisplatin ( Figure 11A ). Based on the targeted genes obtained from the DrugBank database ( Table S7 ), we further found that 15 genes targeted by these drugs had increased expression in the high-risk group ( Figure 11B ). The above evidence indicated that MRMGPS might help clinicians to select personalized drugs and treatments according to different patients.

We further evaluated the affinity of targeted drugs to genes and used MOE for molecular docking to explore the internal relationship between MRMGPS and 10 anti-PCa drugs. The results demonstrated that 51 out of 75 minimum binding energy of drug-gene pairs were less than -6.0 kcal/mol, except that cisplatin could not be docked ( Figures 11C, D ), indicating that the majority of targeted drugs and MRMGPS-based genes maintained extreme stability of the intersectional combination. These findings revealed that there might be a certain relationship between genes and drugs, which provided a certain theoretical value for the study of drugs for PCa.