Integrated analysis of single-cell and bulk RNA sequencing identifies a signature based on macrophage marker genes involved in prostate cancer prognosis and treatment responsiveness

In the tumor microenvironment, tumor-associated macrophages (TAMs) interact with cancer cells and contribute to the progression of solid tumors. Nonetheless, the clinical significance of TAM-related biomarkers in prostate cancer (PCa) is largely unexplored. The present study aimed to construct a macrophage-related signature (MRS) for predicting PCa patient prognosis based on macrophage marker genes. Six cohorts comprising 1056 PCa patients with RNA-Seq and follow-up data were enrolled. Based on macrophage marker genes identified by single-cell RNA-sequencing (scRNA-seq) analysis, univariate analysis, least absolute shrinkage and selection operator (Lasso)-Cox regression, and machine learning procedures were performed to derive a consensus MRS. Receiver operating characteristic curve (ROC), concordance index, and decision curve analyses were used to confirm the predictive capacity of the MRS. The predictive performance of the MRS for recurrence-free survival (RFS) was stable and robust, and the MRS outperformed traditional clinical variables. Furthermore, high-MRS-score patients presented abundant macrophage infiltration and high-expression levels of immune checkpoints (CTLA4, HAVCR2, and CD86). The frequency of mutations was relatively high in the high-MRS-score subgroup. However, the low-MRS-score patients had a better response to immune checkpoint blockade (ICB) and leuprolide-based adjuvant chemotherapy. Notably, abnormal ATF3 expression may be associated with docetaxel and cabazitaxel resistance in PCa cells, T stage, and the Gleason score. In this study, a novel MRS was first developed and validated to accurately predict patient survival outcomes, evaluate immune characteristics, infer therapeutic benefits, and provide an auxiliary tool for personalized therapy.


Introduction
Prostate cancer (PCa) is the world's second most diagnosed cancer type and the leading cause of cancer-related deaths in men (Vietri et al. 2021). Currently, radical prostatectomy (RP) and radiotherapy with or without androgen deprivation therapy (ADT) remain the standard therapeutic regimens for clinically localized PCa (Cornford et al. 2021;Mohler et al. 2019). Although most PCa patients initially respond to treatment, many eventually progress to castration-resistant prostate cancer (CRPC) or neuroendocrine prostate cancer (NEPC) (Ge et al. 2020). The rate of biochemical recurrence (BCR) following RP has been estimated to be 17-33% (Matsumoto et al. 2018;Ward et al. 2003). BCR may predict disease recurrence and subsequent lethal metastases, but BCR often predates other signs of clinical progression by several years (Paller and Antonarakis 2013). Hence, there is an urgent need to determine biomarkers that predict potential recurrence before therapy is initiated. Most studies have focused on cell-autonomous changes in CRPC, with fewer studies on the tumor microenvironment. Among various malignancies, TAMs are the most abundant immune subset, especially in PCa (Bingle et al. 2002;Gollapudi et al. 2013;Hu et al. 2015). The tumor-promoting M2 phenotype is often more commonly observed than the cytotoxic M1 phenotype and is used to reflect the clinical features of advanced PCa (Chanmee et al. 2014). TAMs are known to promote tumor growth and resistance to therapy by releasing cytokines, inhibiting immune surveillance, enhancing angiogenic mediators, and other mechanisms (Izumi and Mizokami 2019;Karan et al. 2009). Mechanistically, macrophage-derived CCL5 can mediate STAT3-dependent epithelial-mesenchymal transition, causing drug resistance and metastasis in PCa (Ma et al. 2021). Research has shown that CSF1R inhibitors may contribute to the therapeutic resistance of PCa as macrophage-targeting agents (El-Kenawi et al. 2021).
Currently, growing scRNA-seq and microarray research associated with tumor heterogeneity and mechanisms has resulted in the emergence of multiple potentially predictive biomarkers for PCa (Baures et al. 2022;Mukherjee and Sudandiradoss 2021;Schnepp et al. 2020). Studies have found that a signature based on B-cell marker genes for lung adenocarcinoma could be used to effectively predict patient survival and immunotherapy (Song et al. 2022). Furthermore, several studies have found that PTGS2 is related to the neurotropic-related genes, and that PTGFR can be used as a potential transcriptome biomarker for predicting the progression of prostate cancer, especially the disease stage (Alkhateeb et al. 2019;Liu et al. 2021a, b). Although the understanding of the role of TAMs in antitumor immunity is increasing, it is still limited. Therefore, it is necessary to use and optimize machine learning algorithms to incorporate the molecular characteristics of immune cells and their correlation with prognosis and immunotherapy prediction into preclinical models, which may provide important biomarkers for biological verification.
Here, a novel MRS model was developed and validated using six independent public datasets to predict patient prognosis and recurrence as well as immune characteristics, drug sensitivity, and immunotherapy response. The relationship between this model and the advantages of leuprolide-based adjuvant chemotherapy (ACT) was analyzed, and drug sensitivity in PCa was assessed. This study may help to explore the biological characteristics and clinical significance of TAMs in cancer, form individualized treatment strategies, and provide a new perspective on cancer immunity.

Curation and preprocessing of publicly available datasets
The genomic expression and matched clinicopathological annotations of prostate cancer samples were retrieved from The Cancer Genome Atlas Prostate Adenocarcinoma (TCGA-PRAD) (https:// portal. gdc. cancer. gov/) and the Gene Expression Omnibus (GEO) database (https:// www. ncbi. nlm. nih. gov/ geo/) (GSE116918, GSE70768, GSE70769, GSE46602 and GSE21032/Memorial Sloan Kettering Cancer Center (MSKCC)). In addition, a scRNA-seq dataset comprising three PRAD samples from GSE153892 was employed to identify macrophage marker genes in different cell clusters. Next, drug-related PCa cell line data were downloaded from GEO datasets (GSE36135-GPL571, GSE33455, and GSE158494) to assess important gene expression in docetaxel-and cabazitaxel-resistant PCa cells. All counts or fragments per kilobase of transcript per million (FPKM) were then converted to transcripts per kilobase million (TPM) for further unbiased comparison. Probe IDs were annotated to gene symbols according to the corresponding GENCODE (Homo sapiens GRCh38). The probes with the same gene were expressed as average values. Finally, a total of 1056 patients with available clinical annotations were included in the subsequent analysis and integrated using the ComBat algorithm, which eliminates potential batch effects across datasets. The baseline information of the collected cohorts is provided in Supplementary Table 1.

Distinguishing and extraction of macrophage marker genes
The scRNA-seq analysis was conducted by using the "Seurat" and "SingleR" packages (Aran et al. 2019;Hao et al. 2021). All other parameters were run with default values. Cells with more than 5% of the mitochondrial gene were removed. Principal component analysis (PCA) dimensionality reduction was performed using 2000 highly variable genes and the optimal clusters were selected for visualization based on the "tsne" package. The criteria used to determine marker genes were |log2FC|> 0.25 and an adjusted P value < 0.05.

Development of the MRS signature
Machine learning algorithms, including Lasso (Gao et al. 2010), stepwise Cox, random survival forest (RSF) (Taylor 2011), generalized boosted regression modeling (GBM), partial least squares regression for Cox (plsRcox), CoxBoost, survival support vector machine (survival-SVM), and supervised principal components (SuperPC) were integrated to build a consensus signature. The "survival" and "survminer" packages were used to fit univariate Cox regression in the meta-cohort, and a P value of 0.05 was considered RFSassociated genes. The Lasso method was performed via the glmnet package, where all crossover mRNAs were penalized to prevent overfitting. The above analysis was adopted to better choose stable hub genes. To make full use of the data, we adopted leave-one-out cross-validation (LOOCV), mainly to prevent overfitting and to evaluate the generalization ability of the model. Consequently, we employed integrative algorithms to establish a predictive model based on the LOOCV framework in the TCGA-PRAD cohort and applied all models to five validation cohorts (GSE116918, GSE70768, GSE70769, GSE46602, and MSKCC). After calculating the concordance index (C-index) across all validation cohorts, the optimal model was selected according to the highest average C-index.

Somatic mutation profile analysis
According to mutation data of the TCGA-PRAD samples, we calculated tumor mutational burden (TMB) using "maftools" to compare the relationship between the TMB and MRS scores. Additionally, OncoPrint plots were generated to display differences in the frequency of alterations between the high/low-MRS-score subgroups. Microsatellite instability (MSI) was inferred using the PreMSIm method.

Evaluation of the immunological characteristics
With transcriptome-based algorithms, ESTIMATE (Yoshihara et al. 2013) was used to compute the immune/stromal scores to estimate the overall immune infiltrations, and single-sample gene set enrichment analysis (ssGSEA) (Hänzelmann et al. 2013) and CIBERSORT (Newman et al. 2015) were used to quantify the relative proportions of immune cells. In addition, 20 immune checkpoint molecules were retrieved from the studies of Auslander et al. (2018) and Hu et al. (2021). Patients' potential response to immunotherapy was inferred by the tumor immune dysfunction and exclusion (TIDE) score.

Functional enrichment analysis
Gene annotation enrichment analysis was performed using the "clusterProfiler" and "fgsea" packages. Differentially expressed genes (DEGs) were identified with a threshold of |log2 FC|> 1 and P < 0.05 by the R package "limma" and clustered into various Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Biological Processes of Gene Ontology (GOBP). Using the file (msigdb. v7.4.symbols.gmt) as a reference, gene set enrichment analysis (GSEA) was performed to analyze various potential signaling pathways. In GSEA, the biological pathways were significantly enriched when the nominal P value was < 0.05 and the FDR < 0.25. Protein-protein interaction (PPI) networks were generated using the "STRINGdb" package to describe the interactions between representative DEGs and visualized using Cytoscape (3.8.2).

Prediction of therapeutic sensitivity
A total of 2191 compounds from three drug sensitivity databases (GDSC2, CTRPv2.0, and PRISM) were accessed through the Cancer Dependency Map (DepMap) portal (https:// depmap. org/ portal/). Lower IC50 and AUC values represent higher sensitivity to compounds. The IC50 (or AUC) of each compound was predicted using the calcPhenotype function of the "pRRophetic" package. The correlation between IC50, AUC values, and MRS scores for prostate cancer was evaluated by Pearson correlation analysis (GDSC: r < − 0.29; CTRP: r < − 0.39; PRISM: r < − 0.35). Boxplots were used to assess the differential drug response of the signature.

Statistical analysis
All statistical analyses were conducted using R 4.2.1 software. Kaplan-Meier (KM) and log-rank tests were applied to evaluate the survival differences. The optimal cutoff value was defined using the "survminer" algorithm. The comparison of the C-index with different models was carried out using the CompareC package. The ROC and area under the ROC curve (AUC) were used to explore exactly the predictive ability, which was plotted using the pROC and timeROC packages. Heatmaps were generated using the R package "ggcor" and "pheatmap". Decision curve analysis (DCA) was used to capture clinical efficacy, while complex models were attempted to be used as decision-making tools (Kerr et al. 2016;Vickers et al. 2008). Using the Pearson analysis was used to evaluate the corresponding coefficients between two continuous variables. The Wilcox test, Kruskal-Wallis test, or Student's t test was used to compare differences between groups. The chi-square test was applied for categorical variables. P values or adjusted P values less than 0.05 were regarded as statistically significant. Figure 1 depicts the complete work framework of this study. The clinicopathological characteristics of the TCGA cohorts are presented in Table 1. T stage, N stage, prostate-specific antigen (PSA) level, and Gleason score showed differences between the low-and high-MRS-score subgroups (P < 0.001). In this study, 1056 PCa patients were included in a large meta-cohort, and the cross-platform batch effect was first removed using ComBat. Following removal, the clustering across platforms was closer than before ( Fig. 2A, B). After data processing and filtering, three PCa samples were extracted from GSE153892, and the gene expression profiles of 11,527 cells were retained for subsequent analysis. The diagram showed the distribution of ten clusters (Fig. 2C, D). Cellular annotations for each cluster were determined by cross-referencing DEGs with canonical marker genes, and cells in cluster 5 were classified as macrophage cells (Fig. 2E). As a result, 447 macrophage marker genes of PCa were identified. Subsequently, 329 genes were obtained from the meta-cohort.

Integrative construction of the MRS model
Based on the 329 macrophage genes from the meta-cohort, univariate analysis was used to screen 104 prognostic genes with P < 0.05 (Supplementary Table 2). Then, Lasso-Cox analysis was performed to acquire 40 more stable genes that were significantly associated with RFS in PCa patients. Finally, these genes were selected for the machine learningbased integration procedure to construct a consensus MRS. In the training dataset, we fitted nine prediction models through the LOOCV framework and then calculated the C-index of each model in the other five validation datasets (Fig. 3A). Interestingly, the average C-index (0.699) was highest in the RSF model and higher in most of the validation datasets (Supplementary Table 3). Meanwhile, the optimal cutoff value was confirmed to be 13.45. PRAD samples were categorized into low-and high-MRS-score subgroups based on this value. The survival curve in the TCGA-PRAD training dataset indicated that high-MRS-score patients exhibited inferior RFS (P < 0.001) (Fig. 3B). Furthermore, similar trends were detected in the five validation datasets and the meta-cohort combining all samples, indicating that the MRS performed well ( Fig. 1G), respectively. The AUC values for the 3-, 4-, and 5-year survival rates were 0.832, 0.874, and 0.816 in the GSE70769 dataset, respectively ( Supplementary Fig. 1E). These results suggested that the MRS model had significant predictive power. The performance of the MRS was further compared with that of other clinicopathological features for predicting prognosis. The distinct accuracy of the MRS was better than that of other variables, including age, T stage, PSA level, Gleason score, TMB, and MSI ( Fig. 4A-F). However, there was no difference in the MSKCC dataset. Next, the performance of the MRS in the TCGA cohort was assessed to further evaluate its clinical utility. The DCA analysis showed that compared to clinical characteristics, the MRS model was associated with more clinical benefits. As described in Fig. 4G, prediction with all or nonpatient scenarios was more beneficial when the decision probability was greater than 0.1 and less than 0.3. The findings also revealed that the combination of the MRS and pathologic parameters might improve the superior clinical usefulness for PCa patients.

Immune landscape of MRS subgroups
We used the GSEA algorithm to identify the various enriched signaling pathways in the meta-cohort to better understand the underlying mechanisms associated with the MRS in PCa. Signaling pathways such as the cell cycle and DNA replication were specifically active in the high-MRSscore subgroup. In contrast, drug metabolism cytochrome p450, the adipocytokine signaling pathway, and peroxisome were significantly enriched in the low-MRS-score subgroup (Supplementary Fig. 2A). Simultaneously, GOBP analysis was performed, and the results suggested that most genes  Fig. 2B). In particular, a thorough search of the TCGA training dataset revealed that 313 DEGs (|Log2FC|> 1 and P < 0.05) were involved in multiple KEGG immune-related pathways, such as the IL-17 signaling pathway, the TNF signaling pathway, the MAPK signaling pathway, human T-cell leukemia virus 1 infection, cytokinecytokine receptor interaction, and the PI3K-Akt signaling pathway. Afterward, the PPI network was generated using the STRING database to reveal interactions among representative DEGs (containing 46 proteins) that regulated the immune response at a specific level (Fig. 5A).
In combination with the enrichment differences observed above, the relationship between the MRS score and immune cells and checkpoints in meta-PCa samples was investigated. The immune, stromal, and ESTIMATE scores were positively correlated with the MRS score, while tumor purity had significant negative dependencies (Fig. 5B, C; Supplementary Fig. 2C, D). Furthermore, CIBERSORT analysis showed that PCa patients with a high MRS score had significantly greater proportions of memory B cells, resting memory CD4 T cells, regulatory T cells (Tregs), gamma delta T cells, activated NK cells, and M1/M2 macrophages but had lower proportions of plasma cells and T follicular helper cells (Fig. 5D). Simultaneously, a correlation heatmap was generated to illustrate the relationships among the MRS score and 28 infiltrating immune cells from ssGSEA. Pearson analysis was performed in all types of immune cells, and 13 pairs exhibited significantly positive correlations, especially between memory B cells and gamma delta T cells, while activated B cells, CD56 bright natural killer cells, eosinophils, immature dendritic cells, mast cells, and natural killer cells were negatively correlated. In addition, the correlation between the expression of 20 genes (15 checkpointrelated genes and 5 immunological activity-related genes) and the MRS score was analyzed. This analysis further suggested that 12 genes were positively correlated with the MRS score and that CEACAM1, CD40, and PVR were negatively correlated (Fig. 5E). Violin plots also suggested that 11 immune cells were highly infiltrated in the high-MRS-score subgroup, which exhibited a positive correlation with the MRS score ( Supplementary Fig. 2E). Moreover, the expression of most checkpoints was upregulated in the high-MRS-score group, which may be associated with the overall upregulation of immune activity. Details of the expression patterns of checkpoints were depicted in Supplementary Fig. 2F. These genes were all positively correlated with the MRS score. Due to the suppression of critical immune checkpoints that may influence cancer immunotherapy with ICB, some representative molecules (CTLA4, HAVCR2, and CD86) were evaluated and found to have a positive correlation with the MRS score. Patients were categorized into four subtypes based on the MRS score and checkpoints, and the impact of their interaction on RFS in PCa patients was determined using the KM curve. The log-rank test results showed that the MRS score could be used to effectively differentiate the prognosis of patients with conflicting levels of CTLA4 (Fig. 5G), HAVCR2 (Fig. 5H), and CD86 (Fig. 5I). Patients with the worst prognosis had low MRS scores and high expression of checkpoints, while patients with high MRS scores and low expression of checkpoints had the highest survival rates among the four subtypes. The Wilcoxon test results also demonstrated that tumors with low MRS scores had significantly lower TIDE scores, implying a better response to ICB therapy (Fig. 5F). Hence, we hypothesized that the MRS score might be an ideal predictor of ICB therapy for PCa patients and correlate with immune activation.

Somatic mutation characteristics were observed in different MRS groups
Increasing evidence suggests that high mutational load and neoantigen overexpression increase the chances of tumor recognition by the immune system ). Therefore, the maftools approach was used to investigate the possible relationship between mutational load and the MRS score. First, patients were classified into two groups (low and high) based on their MRS score. The OncoPrint plots showed that 10 genes were mutated more frequently in the high-MRS-score subgroups (31%, 14%, 10%, 10%, and 10% for TP53, TTN, ABCA13, CSMD3, and FOXA1 respectively) (Fig. 6A). The mutation rates for SPOP, TTN, TP53, FOXA1, and KMT2D in the low-MRS-score subgroup were 12%, 9%, 8%, 5%, and 5%, respectively. Also, the forest plot showed that ABCA13 was highly mutated in the high-MRS-score subgroup (Fig. 6B). Subsequently, a series of evaluation indicators such as TMB, MSI, and immune checkpoint molecule expression were used to represent the response to immunotherapy. The findings suggested that the expression levels of CD86, CTLA4, HAVCR2, PDCD1LG2, and CXCL9 were significantly elevated as the MRS score increased, and the TMB scores were significantly distributed in the high-MRS-score subgroup ( Supplementary Fig. 2F,   Fig. 3 Construction of the MRS model via the machine learningbased procedure and survival curves. A A total of nine prediction models were developed through the LOOCV framework and the C-index was further calculated for each model across all validation datasets. Kaplan-Meier curves showed that patients with lower MRS scores had better RFS than the ones with higher MRS scores in the TCGA-PRAD (B), MSKCC (C), GSE116918 (D), GSE70768 (E), GSE70769 (F), GSE46602 (G), and meta-cohort (H) Fig. 4 Evaluation of the MRS model. A-F The performance of MRS was compared with other clinical and molecular variables in predicting prognosis. G Comparisons of the clinical utility for the clinical variables and combined with MRS using decision curve in the entire TCGA dataset. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; ns: no significance Fig. 6C). The MRS score was higher in MSI-L/MSS than in MSI-H (Fig. 6D). These results demonstrated that the MRS score might influence higher immunogenicity and heterogeneity of PCa.

Potential therapeutic compounds for PCa patients with a high MRS score
Many studies have been searching for potential druggable Fig. 5 Relationship between the pathways and immunological features in prostate cancer. A The diagram depicted the protein interactions of multiple immune-related pathways. The green backward slash lines represent indirectly; the red solid lines represent direct; the gray dotted lines represent the KEGG pathway. Color bars from red to blue represent the fold change in protein levels from increasing to decreasing. The significance of pathways represented by -log (p value) (limma test) was indicated as most significant in dark purple. B, C The box plots presented that both immune score and stromal score are significantly correlated with MRS subtypes. D Comparison of 22 cellular levels inferred by CIBERSORT algorithm between high/low-MRS score subgroups in metacohort. E A correlation heatmap illustrated the relationships of the MRS score with the immune infiltrating cells and inhibitory immune checkpoints. F The differences in TIDE score between two MRS subgroups. KM curve of four subgroups stratified according to the MRS scores and expression of inhibitory immune checkpoints for CTLA4 (G), HAVCR2 (H), and CD86 (I), respectively. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; ns: no significance targets and compounds for treating of PCa patients resistant to standard chemotherapeutic agents. In the TCGA-PRAD cohort, 21 patients were treated with leuprolide-based ACT. The chi-square test results demonstrated that the low-MRSscore subgroup exhibited a higher response rate than the high-MRS-score subgroup (91% vs. 50%) (Fig. 7A), and the MRS score could also be used to significantly discriminate responders from non-responders to leuprolide-based ACT (AUC = 0.722) (Fig. 7B).
Next, three drug response databases (GDSC, CTRP, and PRISM) were searched to identify effective chemotherapeutic agents for patients with high MRS scores and poor prognosis. First, IC50 values were estimated for each sample in the meta-cohort, and then correlations between the values and the MRS score were calculated. After filtering for significant compounds using the criterion of negative r values with P values less than 0.05, the five compounds annotated with the most negative correlation coefficients were displayed in volcano plots (WIKI4, MIM1, vorinostat, GSK2578215A, and WEHI-539), and their estimated IC50 values were lower in the high-MRS-score subgroup (Fig. 7C, D). Moreover, ten compounds were found to be negatively correlated with MRS scores, namely, CD-1530, teniposide, tigecycline, LY-2183240, and leptomycin B in the CTRP database (Fig. 7E) and ixabepilone, triclabendazole, cabazitaxel, paliperidone, and SC-12267 in the PRISM database (Fig. 7G). Furthermore, the results showed that all compounds presented higher AUC values in the low-MRSscore subgroup (Fig. 7F, H). Finally, we obtained a group of drugs from the three datasets whose beneficial therapeutic Fig. 6 Landscape of somatic mutations. A The top 10 mutated genes are in the high-MRS score and low-MRS score subgroups. B The somatic mutations were compared between the high/low-MRS-score subgroups. C The relationship between TMB and MRS scores in PCa patients. D Comparison of the MRS score between MSI-L/MSS and MSI-H groups. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; ns: no significance values were significantly higher in the high-MRS-score subgroup than the low-MRS-score subgroup (p < 0.05, Fig. 7I). Among them, vorinostat, cabazitaxel, and fludarabine have been widely used to treat specific targets of PCa (Carbonetti et al. 2020;Majera et al. 2019).

Association of genes in the MRS model with drug resistance and clinicopathological features
To explore the mechanism of drug resistance, we extracted genes ranked in the top 10 of the RSF in three external datasets based on the PCa cell line and analyzed the expression differences ( Supplementary Fig. 3A-C). Notably, GSE36135 and GSE33455 dataset analysis revealed that ATF3 expression was upregulated in docetaxel-resistant DU-145 cells (Fig. 8A, B). However, the outcome was the opposite in PC-3 cells. Moreover, GSE158494 dataset analysis suggested that ATF3 expression was upregulated in cabazitaxelresistant cells compared to original and docetaxel-resistant cells (Fig. 8C). According to the cutoff threshold of 2.48 derived from the expression of ATF3 in the TCGA cohort, low-expression samples presented significantly poorer RFS than high-expression samples (P < 0.001, Fig. 8D). Similarly, the expression of ATF3 was negatively correlated with the MRS score (Fig. 8E, Supplementary Table 4). In the TCGA cohort, the expression levels of the NAP1L1, IRS2, CYC1, and CKS2 genes were significantly higher in the high-MRS-score subgroup than in the low-MRS-score subgroup. Other genes (such as ATF3, CST3, and CDKN1A) were highly expressed in the low-MRS-score subgroup and were associated with a lower stage (T1-T2) as well as lower Gleason score, predicting a better prognosis for PCa patients. However, no significant difference in survival was observed based on age (Fig. 8F).

Discussion
Macrophages have received significant attention regarding the prognosis of several cancer types (Pittet et al. 2022). Particularly, the presence of TAM infiltration in the tumor microenvironment has been associated with disease progression following ADT, and preclinical studies have indicated that TAMs promote PCa cell proliferation and migration (Cioni et al. 2020;Han et al. 2020;Nonomura et al. 2011). The heterogeneity of PCa is being revisited with the advent of single-cell technologies. Yu et al. discovered FMO2 as a biomarker for macrophage infiltration and prognosis in epithelial ovarian cancer . Additionally, MS4A6A was found to be a new prognostic biomarker produced by macrophages in glioma patients . However, the underlying mechanism of TAM induction remains unclear, and molecular stratification based on predictive biomarkers for TAMs to guide PCa treatment selection has not yet been implemented in clinical practice. Therefore, broadly exploring the prognostic markers for PCa can guide future clinical management.
In the present study, we initially identified macrophage marker genes from PRAD tissues by scRNA-seq analysis. Furthermore, we employed univariate Cox regression analyses and Lasso to screen 40 candidate genes that were highly correlated with the RFS of patients. We further established an integrative method to generate a consensus MRS with the expression profiles of these genes. An advantage of a complex machine learning algorithm is the capacity to develop better statistical models to forecast RFS across all cohorts to enhance classification performance. A total of nine models were fitted to the TCGA-PRAD database through the LOOCV framework. Subsequently, validation results in five independent cohorts obtained from the GEO dataset suggested that RSF was the best model for stratifying PCa patients into two MRS groups. The strength of the model in terms of reproducible and robust performance in predicting RFS was confirmed in the five external cohorts. Furthermore, the AUC value presented satisfying accuracy for molecular subtyping in all cohorts. Previous studies have shown that various clinical predictors are widely used for prognostication and risk assessment of PCa, such as the Gleason score, PSA level, and stage (Albertsen 2020;Truong et al. 2018;Tyekucheva et al. 2017). Other parameters, such as TMB and MSI status, might potentially have an impact on therapy response and prognosis in addition to the clinicopathological features of the patients (Wong and Yu 2021). The C-index suggested that the MRS signature could significantly improve the accuracy of individual relapse risk assessment relative to these factors, which was not observed in Fig. 7 The MRS predicted therapeutic response for PCa patients. A Differences in response rates between the high-MRS score and the low-MRS score subgroup concerning leuprolide-based ACT. B ROC curves of MRS predict the efficacy of leuprolide in TCGA. For the GDSC database, volcano plot (C) and box plot (D) of Pearson's correlations and significance between MRS scores and estimated IC50 values, and blue points being significantly negatively correlated (P < 0.05, r < − 0.29). E, G The AUC values of compounds in CTRP and PRISM were estimated for each PRAD sample and Pearson correlation was performed between MRS scores and estimated AUC values, and the five compounds with the highest negative correlation coefficients were shown on dotted line plots (CTRP: r < − 0.39; PRISM: r < − 0.35). F, H All estimated AUC values for these compounds were significantly lower in the high-MRS-score subgroup. ***P < 0.001; ****P < 0.0001. I Bar plot showed drugs with differences in pharmacological value between high and low MRS scores in the GDSC (blue), CTRP (yellow), and PRISM (pink) databases, where columns represent P values and rows represent drugs. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; ns: no significance respectively. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; ns: no significance. D KM survival curves showed the difference in RFS between the two groups with high and low ATF3 expression. E In the TCGA cohort, ATF3 displayed a significant negative correlation with MRS scores. F The relationship between gene expression and clinical features in the MRS model the MSKCC dataset. Notably, the decision curve, which synthesized the MRS score with clinical characteristics, added a net benefit rate to current clinical values.
Furthermore, these DEGs within subgroups were also correlated with pathways involved in inflammation, such as the IL-17 signaling pathway, the TNF signaling pathway, the MAPK signaling pathway, the PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, peroxisome, and cell cycle. Recent studies have shown that these variables and pathways are essential for metabolic regulation and macrophage polarization (Wang et al. 2019;Zhao et al. 2020). Next, we discovered that PCa patients with high MRS scores were associated with an increased risk of biochemical recurrence and had a higher immune infiltration (such as Tregs, activated CD8 T cells, activated NK cells, and M1 and M2 macrophages). These observations were in agreement with previous studies Mantovani et al. 2017). We observed that higher MRS scores correlated with higher expression of immune checkpoints (CTLA4, HAVCR2, and CD86), which was shown in survival studies to have a mutually beneficial influence on patient prognosis. With the increase in the MRS score, the expression of checkpoints such as CXCL9 and CD80 also increased. This might be the result of the mechanism by which M1-hot TAMs can recruit T cells through CXCL9 expression (Garrido-Martin et al. 2020). According to recent studies, overexpression of HAVCR2 in T cells can lead to dysfunction of PSA-specific CD8 + T cells, further leading to a poor prognosis in PCa patients (Japp et al. 2015;Piao and Jin 2017). Additionally, CD80 binds to CD28 or CTLA4 to activate T-cell costimulation or initiate T-cell co-inhibition, respectively (Chen and Flies 2013). Targeting immune checkpoints has been considered to improve the therapeutic efficacy of immunotherapy modalities in solid tumors (Jafari et al. 2020). This study showed that patients with a low MRS score benefit more from ICB treatment, but those with high TIDE scores may react poorly to immunotherapy. Therefore, the identification of MRS classifiers that predicted the response to immunotherapies was critical to improving their use in treating patients. The mechanism of how the MRS affects RFS and ICB treatment of PCa patients may require more evidence and discussion.
The high-MRS-score subgroup had a higher mutation frequency of TP53 (31%), and the low-MRS-score subgroup tended to have a higher proportion of SPOP mutations (12%). Studies have reported a mutually exclusive relationship between TP53 and SPOP mutations, both of which are independent prognostic markers for CRPC (Mangolini et al. 2022;Stopsack et al. 2020;Zhou et al. 2022). For example, alterations in TP53 were associated with a decreased dependence of the tumor on AR signaling, which was associated with ARSI, whereas the highest levels of AR activity were found in the SPOP mutant subtype, indicating a good response to androgen receptor signaling inhibitor (ARSI) (Liu et al. 2021a, b). Studies have also confirmed that the mutation allele frequencies of PCa patients after treatment mainly occur in TP53, AR, and SPOP mutations ). In addition to modifying CSF1 signaling, the combination of PTEN and p53 loss increased CXCL17 secretion, which further influenced macrophage recruitment and function (Bezzi et al. 2018;Blagih et al. 2020). ABCA13 was highly mutated in the high-MRS-score subgroup. Overexpression of ABCA13 has been reported to be associated with decreased progression-free survival and decreased sensitivity to temozolomide in glioblastoma. (Dréan et al. 2018). Patients with microsatellite instability-high (MSI-H) responded well to the immunological response (Le et al. 2015). This study showed a better immune response in patients with a low MRS score, which is consistent with both the findings from the TIDE approach. However, as confirmed by various studies, the high-MRS-score subgroup had a higher TMB score, suggesting that TMB is unlikely to be a major determinant of the response to immunotherapy for this disease (Barroso-Sousa et al. 2020;Hellmann et al. 2018;Wood et al. 2020). Considering the present and previous findings, it is concluded that the MRS score could be a reliable biomarker to help filter the main population that will benefit from immunotherapy.
Neoadjuvant leuprolide enhances radiation-mediated apoptosis in PCa cells, and its efficacy in reducing tumor volume and prolonging progression-free survival has been well documented (Sánchez et al. 2021). According to these data, patients with a low MRS score tended to have a higher response rate to leuprolide-based ACT. ROC analysis suggested that the MRS provided better accuracy in predicting leuprolide-based ACT. Additionally, potential drug targets and corresponding compounds for PCa patients with a high MRS score were screened based on the established MRS model, and several promising compounds, such as paclitaxel, vorinostat, cabazitaxel, and fludarabine, were selected from three drug response databases. It was discovered that Zn promoted the chemosensitivity of PCa cells to paclitaxel by inhibiting epithelial-mesenchymal transition and inducing apoptosis, which prolonged the life expectancy of PCa patients (Xue et al. 2019). Currently, poly(ADP-ribose) polymerase inhibitors (vorinostat) can be used clinically to activate genotoxic and proteotoxic stress-response pathways in human PCa (Terrisse et al. 2020). Clinical trial data have revealed that cabazitaxel enhanced survival with manageable side effects in patients with metastatic CRPC (Sternberg et al. 2021). Furthermore, increased production of reactive oxygen species by fludarabine phosphate may optimize treatment regimens for patients with N-MYC-overexpressing NEPC tumors (Elhasasna et al. 2022). As a result, the MRS may have practical implications for the development of targeted drugs and combination immunotherapies for PCa patients.
Considering the impact of the tumor microenvironment on drug resistance, a differential analysis of critical genes was performed in external cohorts of PCa cells. The GSE36135 and GSE33455 cohorts showed high expression of ATF3 in docetaxel-resistant DU-145 cells compared to sensitive or original cells; however, the results were opposite in PC-3 cells. Furthermore, GSE158494 data analysis suggested that ATF3 was highly expressed in cabazitaxelresistant cells compared to original and docetaxel-resistant cells. As an inducer of oxidative stress and inflammation, ATF3 may regulate macrophage-associated host defense (Zhao et al. 2021). For example, in response to chemotherapy, wild-type macrophages exhibit pro-oncogenic activity, whereas ATF3 knockout macrophages exhibit anticancer activity (Middleton et al. 2021). ATF3 has been identified as a tumor suppressor for a major subset of PCa with dysfunctional PTEN and has also been shown to inhibit hormoneinduced prostate carcinogenesis in mice (Wang et al. 2016. Our study further showed that ATF3 overexpression was significantly associated with better clinical outcomes and lower MRS scores. Therefore, ATF3 can be used as a key indicator of apoptosis and drug resistance in PCa cells. In addition, the expression of key genes in the MRS model was highly correlated with clinical features, further enhancing its clinical utility. Although external cohorts were validated in this study, single-or multi-center retrospective research data are needed to verify the predictive power of the MRS model, and the prognostic value of patient needs to be thoroughly validated in prospective research. In addition, further investigations are needed to explore the underlying molecular mechanisms between tumor cells and TAM genes. Subsequent clinical trials are essential to further validate the predictive value of potential drug combinations. Ultimately, these efforts are expected to provide a basis for tailored treatment approaches for prostate cancer.

Conclusions
In conclusion, this study developed and validated a novel MRS by integrating machine learning algorithms to effectively predict the RFS of PCa patients. This signature presented superior risk stratification ability and effective prediction of chemotherapeutic drug sensitivity and immunotherapy response, suggesting its promising future utilization.