Integrative multi-omics analysis unveils stemness-associated molecular subtypes in prostate cancer and pan-cancer: prognostic and therapeutic significance
Journal of Translational Medicine volume 21, Article number: 789 (2023)
Prostate cancer (PCA) is the fifth leading cause of cancer-related deaths worldwide, with limited treatment options in the advanced stages. The immunosuppressive tumor microenvironment (TME) of PCA results in lower sensitivity to immunotherapy. Although molecular subtyping is expected to offer important clues for precision treatment of PCA, there is currently a shortage of dependable and effective molecular typing methods available for clinical practice. Therefore, we aim to propose a novel stemness-based classification approach to guide personalized clinical treatments, including immunotherapy.
An integrative multi-omics analysis of PCA was performed to evaluate stemness-level heterogeneities. Unsupervised hierarchical clustering was used to classify PCAs based on stemness signature genes. To make stemness-based patient classification more clinically applicable, a stemness subtype predictor was jointly developed by using four PCA datasets and 76 machine learning algorithms.
We identified stemness signatures of PCA comprising 18 signaling pathways, by which we classified PCA samples into three stemness subtypes via unsupervised hierarchical clustering: low stemness (LS), medium stemness (MS), and high stemness (HS) subtypes. HS patients are sensitive to androgen deprivation therapy, taxanes, and immunotherapy and have the highest stemness, malignancy, tumor mutation load (TMB) levels, worst prognosis, and immunosuppression. LS patients are sensitive to platinum-based chemotherapy but resistant to immunotherapy and have the lowest stemness, malignancy, and TMB levels, best prognosis, and the highest immune infiltration. MS patients represent an intermediate status of stemness, malignancy, and TMB levels with a moderate prognosis. We further demonstrated that these three stemness subtypes are conserved across pan-tumor. Additionally, the 9-gene stemness subtype predictor we developed has a comparable capability to 18 signaling pathways to make tumor diagnosis and to predict tumor recurrence, metastasis, progression, prognosis, and efficacy of different treatments.
The three stemness subtypes we identified have the potential to be a powerful tool for clinical tumor molecular classification in PCA and pan-cancer, and to guide the selection of immunotherapy or other sensitive treatments for tumor patients.
Cancer poses a significant threat to human health and imposes a substantial burden on the global public health. It is the leading cause of death worldwide [1, 2]. Prostate cancer (PCA) is one of the most common malignant tumors in men, and its incidence has increased significantly in recent years [3, 4]. Despite progress in its treatment, such as surgery, radiotherapy, and chemotherapy, the disease remains a challenge, particularly in cases of castrate-resistant prostate cancer [5, 6]. Immunotherapy, which has achieved excellent therapeutic effects in various tumors, is a revolutionary breakthrough in tumor treatment [7, 8]. However, immune checkpoint blockade therapy (ICB) has limited efficacy in unselected patients with PCA, and only a small subgroup of patients may be sensitive to ICB . Therefore, determining patient subtypes that can benefit from immunotherapy is an urgent problem to be solved. Several large-scale clinical trials are currently ongoing, such as the phase III clinical trials of pembrolizumab (KEYNOTE-991) and nivolumab plus docetaxel (Checkmate 7DX), to explore the benefits of immunotherapy with/without other conventional treatments for PCA patients . In recent years, numerous biomarkers that can affect ICB efficiency have been discovered, such as PD-L1, microsatellite instability, tumor mutational burden (TMB), and TCR polymorphism . However, they are far from being ideal biomarkers for PCA .
The molecular subtyping of tumors can be used to guide the precision diagnosis and treatment of tumor patients. Cellular and molecular characteristics of tumors have shown potential in precision therapy of PCA, such as cancer-associated fibroblasts (CAFs) and their signature genes [12, 13]. Serum prostate-specific antigen (PSA) is the most important marker for PCA; however, it has been criticized for its poor specificity . Existing molecular subtyping methods for PCA were designed for specific clinical applications. For example, the PAM50 method is used to guide androgen deprivation therapy ; the Decipher method is used to guide radiotherapy and surgical treatment . In addition, existing biomarkers have limited clinical utility and are not suitable for guiding ICB treatment [17,18,19]. There is still lacking subtyping methods for PCA, which can effectively and systematically characterize patients from various clinical points of view, e.g., diagnosis, prognosis, recurrence risk, metastasis risk, progression risk, and efficacies of different treatments.
Stemness refers to the self-renewal and differentiation potential of cells. In almost all human malignant tumors, there is a rare subset of cancer cells with stem-like properties, called “cancer stem/stem-like cells” (CSCs) . CSCs are considered the origin cells of tumors and play an important role in the recurrence, metastasis, and treatment resistance of many tumors, including PCA [21,22,23]. They also affect the effectiveness of immunotherapy . Therefore, stemness based subtyping holds potential in personalized management of PCA patients.
In this study, we utilized single-cell RNA-seq (scRNA-seq), bulk RNA-seq, methylation array, and whole exon sequencing datasets to systematically assess stemness differences among PCA patients. By integrating scRNA-seq and bulk RNA-seq, we developed a stemness-based subtyping model consisting of 18 stemness related gene-sets that separated PCA samples into three subtypes with distinctive clinical and molecular characteristics, functional annotations, prognoses, and treatment responses, especially immunotherapy. Furthermore, a subtype predictor including 9 stemness related genes was constructed which showed great performance in tumor diagnosis, predicting ICB and androgen deprivation therapy (ADT) response, metastasis, recurrence, progression, and prognosis.
Stemness scores are closely correlated with clinical and molecular features
Single-cell and bulk RNA-seq analyses show a positive correlation between stemness and PCA malignancy
The workflow of this study is depicted in Additional file 1: Fig. S1 and Methods S1, and all the datasets used in this study are presented in Additional file 1: Table S1. Prostate epithelial cells were extracted from five PCA scRNA-seq datasets [25,26,27,28,29] (Additional file 1: Fig. S2) and the cytoTRACE algorithm  was used to calculate their stemness levels. Consistently, the results showed significantly higher cytoTRACE scores in the malignant epithelium than in the para-cancerous or benign prostate epithelium (Fig. 1a, Additional file 1: Figs. S3a–c). Moreover, among malignant epithelial cells, high-grade PCA (Gleason score [GS] > 7) showed significantly higher stemness scores than low-grade PCA (GS ≤ 7) (Fig. 1b, Additional file 1: Fig. S3d).
In parallel, there are stemness indices calculated by the one-class logistic regression (OCLR) algorithm based on specific stemness probes, including mRNAsi obtained from bulk RNA-seq data and mDNAsi, EREG_mDNAsi, and DMPsi obtained from DNA methylation data . We computed the stemness indices of 553 patients from The Cancer Genome Atlas-Prostate Adenocarcinoma (TCGA-PRAD). These results further confirmed the findings obtained from the scRNA-seq (Fig. 1c, Additional file 1: Fig. S3e). Likewise, these results were consistently replicated in six other PCA bulk RNA-seq datasets [32,33,34,35,36] (Additional file 1: Fig. S3f).
Stemness scores are significantly associated with PCA clinicopathological and molecular features
We divided the patients into two groups based on their stemness indices. Kaplan–Meier (K–M) analysis showed that the high stemness indices group had remarkably poorer overall survival (OS), progression-free interval (PFI), disease-free interval (DFI), and disease-specific survival (DSS) (Fig. 1d, Additional file 1: Fig. S3g–i). Furthermore, we demonstrated the association between stemness and worse prognosis on another six PCA datasets [32, 33, 37,38,39,40] (Additional file 1: Fig. S3j). Based on these results, we conclude that stemness is a risk factor for the prognosis of PCA.
Additionally, we found that mRNAsi is closely related with serum PSA, androgen receptor (AR), tumor purity, pathological T stage (pT), as well as somatic copy number alterations (SCNAs) and TMB (Additional file 1: Fig. S3k, l, Remark A). Furthermore, we observed a negative correlation between stemness and immune infiltration levels (Additional file 1: Fig. S4, Remark B).
Identification of three PCA subtypes based on stemness and their association with prognosis
Development of 18 stemness signatures for subtyping
Since there was a significant correlation between stemness scores and PCA prognosis, we postulated that certain signatures indicating the tumor’s stemness status could be utilized for patient classification. Therefore, we conducted a joint analysis of five scRNA-seq [25,26,27,28,29] (75,350 cells) and seven bulk RNA-seq datasets [32,33,34,35,36] (1641 samples) to develop these signatures. The development process was illustrated in Additional file 1: Fig. S5a. We identified 288 genes from the scRNA-seq data that showed a significant positive correlation with cytoTRACE scores and were significantly overexpressed in both tumors (compared to normal samples) and high-grade PCAs (compared to low-grade PCAs, Fig. 2a, left). Additionally, we identified 220 genes from the bulk RNA-seq data that exhibited a significant positive correlation with mRNAsi and were highly expressed in the tumor samples (Fig. 2a, right). After applying union, gene set enrichment analysis (GSEA), and univariate COX analysis, we identified 18 gene-sets representing risk factors (Additional file 1: Fig. S5b) for subsequent stemness subtyping, including 160 genes. Among these 18 gene sets, 7 were metabolic-related, 8 were cell cycle-related, and the other 3 were MYC targets, mTORC1 signaling, and cellular responses to stimuli, all of which are classic oncogenic signaling pathways [41, 42] (Additional file 1: Fig.S5c, Additional file 2: Data S1).
Identification of three stemness subtypes
Subsequently, based on the ssGSEA scores of these 18 signatures, we used unsupervised hierarchical clustering to classify the 553 samples from TCGA-PRAD into three subtypes (Fig. 2b, c, Additional file 1: Fig. S6a). We defined the group with the highest scores, which contained 56 samples (10.1%), as the “High Stemness” (HS) subtype, while the group with the lowest scores, consisting of 261 samples (47.2%), was named the “Low Stemness” (LS) subtype. The medium score group, comprising 236 samples (42.7%), was referred to as the “Medium Stemness” (MS) subtype, indicating that they were in a transitional state with the potential to transform into either a high or low stemness status. To validate our classification results, we compared the four aforementioned stemness indices across the three subtypes. Compared to LS, mRNAsi, mDNAsi, EREG_mDNAsi, DMPsi scores, and tumor purity were all significantly higher in HS and intermediate in MS (Fig. 2d, Additional file 1: Fig. S6b).
Significant differences in prognosis among patients with three stemness subtypes of PCA
K–M analysis revealed that HS patients had the worst PFI, whereas LS patients had the longest PFI (Fig. 2e, p = 0.00015), indicating that HS patients may experience disease progression earlier after treatment. However, there were no significant differences in OS, DFI, or DSS among the three subtypes (Additional file 1: Fig. S6c). Univariate COX analysis showed that LS and HS were protective factors and risk factors for PFI, respectively (Fig. 2f). Moreover, the multivariate COX analysis further demonstrated that LS was an independent protective factor for PFI (Fig. 2g).
Finally, we further validated the stability and universality of the stemness classification method by using two additional algorithms, non-negative matrix factorization (NMF)  and CensusClusterPlus , along with seven independent PCA bulk RNA-seq datasets [32, 33, 37,38,39,40, 45] (Additional file 1: Fig. S6d–u, Remark C).
Three stemness subtypes have distinct clinical features, mutational events, and functional annotations
Significant differences in certain clinicopathological features
Next, we compared the demographic and clinicopathological characteristics of the three stemness subtypes in patients with PCA. We found that LS had the highest proportion of para-cancerous samples (Fig. 3a, b; LS:MS:HS = 18%:1%:2%, p = 6.6e−11). The proportion of high-grade PCAs (GS > 7), pT3 and pT4 PCAs, and lymph node metastasis (N1) patients gradually increased in the LS, MS, and HS subtypes (Fig. 3a and b). There were no significant differences in body weights among the three subtypes (Fig. 3c). Compared with LS and MS, HS patients were older at diagnosis (Fig. 3c), suggesting that PCA patients diagnosed at an older age may have a relatively higher stemness status. Moreover, Gleason score, PSA, and androgen receptor were all significantly higher in HS patients (Fig. 3c), indicating that HS patients with the highest malignancy may benefit more from ADT treatment.
Significant differences in mutation events
Genome changes are research hotspots; therefore, we performed somatic mutation analysis. TMB, fragment genomic alteration (FGA), and amplification were all significantly higher in the HS subtype (Fig. 3d), suggesting that the HS subtype may be more likely to respond to immunotherapy . The deletion was significantly higher in the MS subtype than in the LS subtype (Fig. 3d). Exon imbalance scores did not differ between stemness subtypes (Fig. 3d). Oncoplot revealed that the HS subtype had the highest rate of genomic alterations (Fig. 3e–g, LS:MS:HS = 66.82%:79.83%:92.73%). Although the top mutated genes of the three subtypes are similar, there are great differences in their mutation rates (Fig. 3e–g). We also observed that TP53 and ABCA13 mutations co-occurred, as did TTN and HERC2 mutations in the HS subgroup (Additional file 1: Fig. S7a). Silencing of TP53, a tumor suppressor gene, and inactivation of ABCA13, a cholesterol transporter protein, may be possible mechanisms contributing to the highly malignant phenotype of HS [47, 48]. In the LS subtype, we found a co-mutation of TTN and FAT3 (Additional file 1: Fig. S7b), while SPOP and TP53 mutations were mutually exclusive in the MS subtype (Additional file 1: Fig. S7c). The proportion of more SCNAs cluster gradually increased in the stemness subtypes (Fig. 3h, LS:MS:HS = 16%:38%:72%, p = 2.7e−14), while the high methylation cluster decreased successively (Fig. 3h, p = 2.6e−8). We also investigated the status of mutations and copy number alterations (CNAs) in common biomarkers of PCA. We found that the proportion of TP53 mutations and CNAs (heterozygously deleted [hetloss] + homozygously deleted [homdel]) increased progressively in the three subtypes (Fig. 3h, p = 0.06 and 0.0098, respectively). Compared with LS, the proportion of PTEN mutations was significantly higher in MS and HS (p = 0.03), and the proportion of CNAs increased in these three subtypes (Fig. 3h, p = 0.0022). Similarly, both the mutations and CNAs of CDKN1B and CNAs of RAD51C, FAM175A, CHD1, RB1, FANCC, and SPOPL were significantly higher in HS than in LS (Additional file 1: Fig. S7d). However, we found no significant differences in the mutations and CNAs of BRCA1 and BRCA2 between the three subtypes (Additional file 1: Fig. S7d). Understanding the mutation landscape of the aforementioned stemness subtypes is beneficial for uncovering potential mechanisms underlying tumor development, and provides a basis for discovering potential therapeutic targets and biomarkers.
HS enriches oncogenic signaling pathways
Subsequently, we performed gene set variation analysis  (GSVA), GSEA , and ingenuity pathway analysis (IPA) to investigate the functional annotations, signaling pathways, and underlying mechanisms associated with the PCA stemness subtypes. As shown in Additional file 1: Fig. S8, the LS subtype mainly enriched nonspecific pathways. In contrast, HS exhibited enrichment in cell cycle-related signaling pathways , but IPA showed that the activity of these signaling pathways was suppressed. MS-enriched signaling pathways were similar to those of HS, but with slightly lower expression levels. See Additional file 1: Remark D for more details.
Three stemness subtypes have different treatment sensitivity and TIME patterns
Three stemness subtypes retain sensitivity to specific drugs
Given that the three stemness subtypes have unique functional pathways, we used the oncoPredict package  to predict the sensitivity of the three stemness subtypes to drugs and quantified sensitivity using half-maximal inhibitory concentrations (IC50). For the selection of conventional drugs, HS and MS subtypes were more sensitive to ADT (bicalutamide, abiraterone) (Fig. 4a, Additional file 1: Fig. S9a); for chemotherapy drugs, HS was more sensitive to taxanes (paclitaxel, docetaxel) (Fig. 4b, Additional file 1: Fig. S9b), etoposide (Fig. 4c), and gemcitabine (Fig. 4d), but resistant to mitoxantrone (Fig. 4e) and platinum drugs (Fig. 4f). In terms of recommended drug selection, HS demonstrated greater sensitivity to epidermal growth factor receptor (EGFR) inhibitors such as sunitinib (Fig. 4g), sorafenib, and imatinib (Additional file 1: Fig. S9c, d), while showing lower sensitivity to cabozantinib, afatinib, and erlotinib (Additional file 1: Fig. S9e–g). Additionally, HS showed increased sensitivity to poly (ADP-ribose) polymerase (PARP) inhibitors, including olaparib and talazoparib (Fig. 4h, Additional file 1: Fig. S9h), as well as to BI-2536 (Fig. S9j), histone deacetylase (HDAC) inhibitors (tacedinaline) (Additional file 1: Fig. S9k), and TGF-β receptor inhibitors (SB-431542) (Additional file 1: Fig. S9l). Conversely, MS exhibited greater sensitivity to anti-EGFR (cetuximab) (Additional file 1: Fig. S9m) and LS was more responsive to cyclin-dependent kinase (CDK) inhibitors (AZD5438) (Additional file 1: Fig. S9n).
Three stemness subtypes have distinct TIME patterns and immunotherapy responsiveness
As demonstrated above, the stemness scores exhibited an inverse correlation with immune infiltration in PCA (Additional file 1: Fig. S4). To further investigate the relationship between stemness subtypes and their response to immunotherapy, we examined the TIME patterns. Our results indicated that the LS subtype had the highest Stromal, Immune, and ESTIMATE scores compared to the MS and HS subtypes (Fig. 4i). This indicated that the LS subtype had the most abundant immune infiltration. Moreover, we observed a significant association between the LS subtype and the High.immu subtype, whereas the MS subtype was significantly associated with the Low.immu subtype (Fig. 4j, Additional file 1: Fig. S9o).
We utilized the CIBERSORT algorithm  to assess immunocyte infiltration levels in the 553 PCA samples. Our findings revealed that CD4T, M1 and M2 macrophages had the highest abundance, while mast and plasma cells showed the lowest infiltration in HS (Fig. 4k). Moreover, Spearman correlation analysis showed a dramatically negative correlation between the stemness indices and the majority of immune checkpoint molecules (Fig. 4l), and most molecules were least expressed in HS (Additional file 1: Fig. S9p). These differences in TIME patterns and immune checkpoint molecule expression could potentially affect the efficacy of immunotherapy. Therefore, we applied the TIDE algorithm  to estimate the response of the patients to immunotherapy. The results showed that the proportion of responders in the HS was higher than that in the LS and MS subtypes (Fig. 4m, Additional file 1: Fig. S9q, p = 0.08). Additionally, the submap analysis  revealed that HS was the most sensitive to anti-PD-L1 (Fig. 4n).
Construction and validation of stemness subtype predictor
Construction of stemness subtype predictor
To facilitate the clinical application of our findings, we developed a stemness subtype predictor with high sensitivity and specificity, using the process outlined in Additional file 1: Fig. S10a (see Additional file 1: Methods S1 for details). We performed weighted gene co-expression network analysis (WGCNA)  (Fig. 5a, Additional file 1: Fig. S10b–g), protein–protein interaction (PPI) analysis (Fig. 5b), and Venn diagram plotting in sequence, and ultimately obtained 9 most critical stemness marker genes (Fig. 5c, Additional file 1: Table S2). We observed that these genes were significantly overexpressed in tumors (compared to normal, Additional file 1: Fig. S10h) and high-grade PCAs (compared to low-grade, Additional file 1: Fig. S10i).
We conducted a literature review on the research status of genes in PCA. CDK1 , KIF4A , TPX2 , BUB1 , and TOP2A  have been reported to promote PCA progression. However, there is a lack of evidence regarding the roles of SKA3, DLGAP5, NCAPG, and HMMR in PCAs. We collected 60 PCA samples and performed immunohistochemistry (IHC) to verify the expression of these four proteins. Our findings showed that their expression gradually increased in benign prostate tissues, and in low-grade and high-grade PCAs (Fig. 5d).
Predictor can effectively distinguish stemness subtypes
Consistently, 100 MLs and unsupervised hierarchical clustering based on the 9-gene stemness subtype predictor showed that our predictor is highly effective in accurately distinguishing the stemness subtypes and can be used as an alternative clinically practical marker panel instead of the 18 stemness signatures (Additional file 1: Fig. S11a–d, Remark E).
Stemness subtype predictor has excellent performance in predicting malignancy, recurrence, progression, metastasis, and treatment response
Subsequently, we investigated the performance of our predictor in distinguishing malignancy and predicting tumor recurrence, progression, metastasis, and ADT efficacy. In TCGA-PRAD, GSE21034 , and GSE70770  datasets, the average area under the receiver operating characteristic curves (AUCs) of the top 10 MLs for identifying benign and malignant tumors were 89.8%, 89.5%, and 89.4%, respectively (Additional file 1: Fig. S11e). For the prediction of ADT response, the top 10 MLs had mean AUCs of 88.6%, 80.5%, and 76.5% in GSE6811 , GSE28680 , and Alumkal_2020 , respectively (Fig. 5e, Additional file 1: Fig. S11f). For the prediction of biochemical recurrence, the mean AUCs of the top 10 MLs in prostate_dkfz_2018 , GSE70770 , and GSE21034  were 78.9%, 74%, and 72.8%, respectively (Fig. 5e, Additional file 1: Fig. S11f). In predicting metastasis, the mean AUCs were 89.4% in GSE21034 , 74.3% in GSE116918  (Fig. 5e), and 66.3% in GSE46691  (Additional file 1: Fig. S11f). The average AUC for the prediction of progression was 65.8% (Fig. 5e). Detailed data of the above results can be found in Additional file 2: Data S2.
Stemness subtype predictor can serve as a prognostic risk model
We calculated stemness-related risk scores (SRS) for each patient based on the above nine genes. Samples were then categorized into high and low SRS groups, and K–M analysis showed that the high SRS group had significantly worse PFI, OS and DFI than the low SRS group (Fig. 5f, Additional file 1: Fig. S11g–i). Univariate COX analysis revealed that SRS, Gleason score, pT, pN, and purity were the risk factors for PFI (Additional file 1: Figure S11j). Multivariate COX analysis further confirmed that SRS was the most significant independent prognostic indicator of PCA (Fig. 5g, HR = 2.04, p = 0.005). Furthermore, receiver operating characteristic (ROC) analysis demonstrated that SRS can serve as an effective clinical prognostic biomarker (Additional file 1: Fig. S11k–o). The K–M analysis based on another six independent PCA datasets [32, 37, 39, 67,68,69] further validated the above results (Additional file 1: Fig. S11p–u, Remark E).
Conserved stemness subtypes in pan-tumors
Three stemness subtypes with significantly different prognoses are prevalent in pan-tumors
To determine whether our stemness subtypes were conserved in pan-tumors, we downloaded four pan-tumor datasets: TcgaTargetGtex, PCAWG, ICGC, and GSE2109. Using the 18 stemness signatures, we consistently classified the datasets into three stemness subtypes, each with significantly different prognoses (Additional file 1: Fig. S12a–h). This suggests that stemness subtypes have a general commonality in pan-tumors.
Furthermore, we confirmed that the stemness subtype predictor can effectively distinguish pan-tumor stemness subtypes, and SRS can predict the clinical prognosis of pan-tumor patients (Additional file 1: Fig. S12i–s, Remark F; Additional file 2: Data S2).
HS patients are more inclined to benefit from immunotherapy
Here, we investigated the correlation between stemness subtypes and ICB responsiveness in pan-tumors. We collected RNA-seq data from 2641 pretreatment samples from patients across 36 datasets and 12 tumor types who received ICB therapy [70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102]. Figure 6a clearly showed that these samples could be classified into three subtypes based on 18 stemness signatures. Hypergeometric test showed that stemness subtypes were significantly correlated with ICB responsiveness (Fig. 6b). The response rate increased progressively across the three subtypes (Fig. 6c). These results were further corroborated by K–M analysis, submap analysis , and separate clustering analysis of each cancer type (sample size > 50) (Additional file 1: Fig. S13a–g, Remark G).
Subsequently, univariate logistic regression analysis was performed based on the RNA-seq data of pre-ICB treatment samples to evaluate the effect of the three stemness subtypes on immunotherapy outcomes. Our findings revealed LS as a risk factor and HS as a promoting factor for response (Fig. 6d). Similar observations were also made in separate analyses of bladder cancer (BLCA), breast cancer (BRCA), non-small cell lung cancer (NSCLC), melanoma (SKCM), and head and neck squamous cell carcinoma (HNSCC) (Additional file 1: Fig. S13h). Nonetheless, the stemness subtypes did not significantly affect immunotherapy in clear cell renal cell carcinoma (ccRCC) (Additional file 1: Fig. S13h). We further compared the stemness subtypes with other ICB predictive biomarkers using the IMvigor210 dataset. The results showed that the stemness subtype was an effective predictor of ICB responsiveness (Additional file 1: Fig. S13i–l, Remark G).
Next, we attempted to further explore the effect of immunotherapy on stemness. We collected bulk RNA-seq data of paired samples from the same patients before and on/post ICB treatment [80, 82, 86, 93, 96, 98,99,100, 102]. Using hierarchical clustering to distinguish these samples into three subtypes (Additional file 1: Fig. S14a), we found that the proportion of HS decreased significantly after ICB therapy, while the proportion of LS increased significantly (Additional file 1: Fig. S14b). Further distinguishing between responders and non-responders, we found that most of the HS patients who responded to ICB treatment shifted to LS after ICB treatment (Fig. 6e, Additional file 1: Fig. S14c), while most of the HS patients who were resistant remained unchanged after treatment (Fig. 6f, Additional file 1: Fig. S14d). These observations were consistent across analyses of HNSCC, esophageal adenocarcinoma (EAC), and SKCM (Additional file 1: Fig. S14e–g). We also validated our findings using a single-cell dataset of HNSCC, which demonstrated a significant decrease in cytoTRACE scores following ICB treatment (Fig. 6g).
We applied single-cell and bulk RNA/DNA-seq technologies to evaluate the stemness status of patients based on stemness signatures. This allowed us to subtype the patients into three stemness subtypes. Specifically, HS patients are sensitive to androgen deprivation therapy, taxanes, and immunotherapy and have the highest stemness, malignancy, TMB levels, worst prognosis, and immunosuppression. LS patients are sensitive to platinum-based chemotherapy but resistant to immunotherapy and have the lowest stemness, malignancy, and TMB levels, best prognosis, and the highest immune infiltration. MS patients represent an intermediate status of stemness, malignancy, and TMB levels with a moderate prognosis. Our 9-gene stemness subtype predictor demonstrated high sensitivity and specificity, and could be easily used to identify stemness subtypes of patients through real-time quantitative PCR or IHC. This simple procedure provides a powerful tool for clinical tumor typing and treatment selection.
PCSCs have been shown to play a significant role in the occurrence, development, treatment resistance, and recurrence of PCA . In our study, the HS subgroup exhibited a significantly worse prognosis (Fig. 2e), significantly higher Gleason score, and expressed the highest level of serum PSA (Fig. 3c), indicating the highest malignancy in HS patients. It is essential to study the signaling pathways and functional annotations associated with PCSCs to understand the molecular mechanisms of carcinogenesis and to identify potential drug targets. Several signaling pathways are closely related to the regulation of PCSCs, such as the PI3K/AKT/mTOR pathway and the Wnt/β-catenin pathway [103, 104]. Studies have shown that PCSCs have a relatively slow proliferation rate compared to ordinary PCA cells and exist in a quiescent status [105, 106]. GSVA and GSEA revealed that cell cycle-related signaling pathways, including cell cycle checkpoints, were enriched in HS (Additional file 1: Fig. S8a, b). However, IPA showed that the cell growth, proliferation, and development of HS were significantly inhibited (Additional file 1: Fig. S8e). These results suggested that while the signaling pathways and molecules related to the cell cycle are upregulated in HS, their proliferative activity is suppressed or static. This phenomenon may be attributed to the self-renewal and multidirectional differentiation ability of PCSCs, such as negative feedback regulation after excessive proliferation (enrichment of cell cycle checkpoints) and activation disorders of cell cycle-related proteins [106,107,108,109].
The PCA treatment continues to be a major research focus. Radical prostatectomy and radiation therapy, with or without ADT, are the standard treatments for localized PCA. ADT remains the primary treatment for advanced disease [10, 110]. Although many treatments initially eradicate cancer cells, cancer often recurs owing to the presence of resistant PCSCs . Therefore, targeting CSCs is critical for preventing tumor recurrence. The significance of molecular subtyping in tumors lies in its ability to help doctors better understand the biological characteristics, molecular mechanisms, development trends, and treatment responses of tumors, and to provide personalized treatment plans for patients. Existing molecular subtyping methods for PCA include the PAM50 method , which is used to guide ADT, and the Decipher method , which is used to guide radiotherapy and surgical treatment. By contrast, our stemness subtyping method can be used to guide various treatments. We found that different stemness subtypes exhibited different sensitivities to drugs. HS was more sensitive to ADT (bicalutamide), PARP inhibitors (olaparib), EGFR inhibitors (sunitinib and sorafenib), immunotherapy, and chemotherapeutic drugs such as docetaxel, etoposide, and gemcitabine. In contrast, LS was more sensitive to platinum drugs, erlotinib, and CDK inhibitors (AZD5438) (Fig. 4, Additional file 1: Fig. S9). Previously, the effectiveness of immunotherapy for PCA has been limited. However, with the development of immunology and cutting-edge molecular diagnostic tools, immunotherapy is expected to become a viable treatment option for PCA . Unlike tumors, such as melanoma, bladder cancer, and NSCLC, which are highly responsive to immunotherapy and characterized by infiltrating lymphocyte proliferation, PCA is considered a “cold” tumor with an immunosuppressive TME [9, 112]. Although the effectiveness of ICB monotherapy for PCA is limited, combined strategies with other standard treatments (ADT, chemotherapy, PARP inhibitors, radium-223, and tyrosine kinase inhibitors) have shown some positive effects . Overall, these findings underscore the importance of molecular subtyping for guiding cancer treatment. By tailoring therapies to the specific molecular characteristics of tumors, doctors can improve treatment outcomes and help patients achieve the best possible outcomes.
The significance of pan-cancer research lies in the application of diagnosis and treatment to more tumors through cross-tumor similarities . We validated the conserved characteristics of the three subtypes in pan-tumors and their responsiveness to immunotherapy using 28,381 pan-tumor samples and 2641 ICB pretreatment samples. This indicates that our findings have implications beyond PCA, and can potentially benefit a wider range of patients. Furthermore, developing a predictive model using multiple MLs and datasets can enhance the model’s generalization ability and prevent overfitting or underfitting issues, which is an effective way to improve the accuracy and robustness of the predictors. We jointly developed a 9-gene stemness subtype predictor with high sensitivity, specificity, and excellent generalization ability using four datasets and 76 machine learning algorithms. Significantly, this predictor can be further developed into a kit for clinical application. We believe that this predictor has great potential for clinical application, as it offers rapid and reliable molecular diagnosis and prognosis for patients with PCA and guides personalized treatment decisions. Moreover, this predictor can facilitate the enrollment of PCA patients into clinical trials for immunotherapy or other targeted therapies based on their stemness subtype, and it may also be applicable to other cancer types with the same stemness subtypes as PCA.
Although this study yielded valuable insights, it is important to acknowledge some of its limitations. First, the sample size for our real-world validation was only 60, and we did not include any pan-cancer samples. Expanding the sample size and including pan-cancer samples for verification in future studies are essential. Second, there is a lack of RNA-seq data of PCA treated with ICB. Although the pan-tumor ICB treatment cohort validated the relationship between stemness subtypes and immunotherapy responsiveness, further validation is still needed for PCA immunotherapy datasets.
In conclusion, tumor molecular typing is of great significance for understanding how cancer develops and progresses, as well as for guiding clinical treatment and the development of new anticancer drugs. By categorizing PCA patients into three stemness subtypes, we can systematically characterize patients from various points of view, e.g., stemness, prognosis, clinical pathological features, mutation patterns, malignancy degree, immune infiltration levels, and efficiencies of different treatments, including immunotherapy. This classification method is also applicable to pan-tumor analyses. Furthermore, the 9-gene stemness subtype predictor we developed is expected to be a clinically useful tool for precision oncology. Significantly, our method provides a pipeline for the development of cancer classification that can be applied to various tumors based on different research hotspots.
Patients and datasets collection
Five PCA scRNA-seq datasets [25,26,27,28,29], 20 PCA bulk RNA-seq datasets [32,33,34,35,36,37,38,39,40, 45, 62,63,64,65,66,67,68,69], four bulk-RNA-seq datasets of pan-tumors, one ICB-treated bulk RNA-seq cohort of pan-tumors (consisting of 36 independent datasets) [70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102], a pluripotent stem cell (PSC) expression matrix, and DNA methylation, somatic mutation, and copy number alteration data from TCGA-PRAD were used in this study. The database used in our study included the MSigDB  (v2023.1.Hs, http://www.gsea-msigdb.org/gsea/index.jsp), and STRING  (v11.5, https://cn.string-db.org/). See Additional file 1: Methods S1 and Table S1 for more details.
After obtaining patient consent and approval from the institutional research ethics committee, we collected paraffin-embedded tissue sections of prostate cancer and benign prostatic hyperplasia (BPH) from the Pathology Department of the Shanghai Sixth People’s Hospital.
The cytoTRACE (v0.3.3) package  was used to calculate cytoTRACE scores for stemness evaluation in single-cell RNA-seq data. The OCLR machine learning algorithm  quantified stemness levels using mRNAsi, mDNAsi, EREG_mDNAsi, and DMPsi scores in bulk RNA-seq data. Higher scores indicate greater stemness levels. In K–M analysis, the surv_cutpoint function in the Survminer (v0.4.9) package calculated the optimal cutoff point to divide patients into high and low stemness groups. See Additional file 1: Methods S1 for further details.
Unsupervised hierarchical clustering was used for clustering analysis based on stemness signatures, implemented using the hclust function. Consensus clustering and NMF were employed using the ConsensusClusterPlus  (v1.58.0) and NMF  (v0.24.0) packages. The optimal number of clusters was determined using the consensus heatmap, cumulative distribution function (CDF) curves, and the proportion of ambiguous clustering algorithm (PAC)  (Additional file 1: Fig. S6a). See Additional file 1: Methods S1 for further details.
Survival analysis was conducted using the Survival (v3.4.0) and Survminer (v0.4.9) packages. The Kaplan–Meier method plotted the survival curve, and the log-rank test compared survival differences between groups. The Cox proportional hazards model investigated covariate effects on survival time, calculating the risk ratio and confidence interval.
Additionally, the timeROC package (v0.4) was used to estimate the time-dependent ROC curve and AUC, which allowed us to evaluate the prognostic ability of SRS.
Identification of stemness marker genes and developing of signatures for stemness classifications
The identification process is illustrated in Additional file 1: Fig. S5a. We utilized scRNA-seq data to identify stemness marker genes in malignant and high-grade PCA epithelial cells, and used bulk RNA-seq and DNA methylation data to identify stemness marker genes in malignant PCA samples. Our approach involved analyzing multiple datasets, performing Spearman correlation analysis and DE analysis, and selecting genes that were positively correlated with stemness scores and were upregulated in malignant and high-grade PCA. We then analyzed gene sets on MSigDB  with the stemness marker genes obtained and identified significant risk signaling pathways using univariate COX analysis for subsequent stemness classification. See Additional file 1: Methods S1 for further details.
Paraffin sections underwent dewaxing, antigen retrieval, and serum blocking. They were incubated with primary antibodies overnight at 4 °C, secondary antibodies and SABC for 30 min at 37 °C. Sections were stained with DAB and counterstained with hematoxylin. The primary antibodies used were the NCAPG Polyclonal antibody (24563-1-AP, Proteintech, Wuhan, China), HMMR-specific polyclonal antibody (15820-1-AP, Proteintech, Wuhan, China), HURP polyclonal antibody (12038-1-AP, Abcam, Proteintech, Wuhan, China), and SKA3 antibody (SC-390326, Santa, California, USA).
Numerical variables were compared using t-tests or ANOVA, categorical variables using χ2, Fisher’s exact or Kruskal–Wallis tests. Non-normally distributed variables were compared using non-parametric tests. Correlations were evaluated using Pearson or Spearman tests. Survival differences were compared using the log-rank test. Confidence intervals (CIs) were reported as 95% and significance was set at P < 0.05. Analyses were performed using R (v4.2.1), Python (v3.10), and Origin 2022 software.
More methods and details can be seen in Additional file 1: Methods S1.
Availability of data and materials
The original contributions presented in the study are included in the article/Additional files. Further inquiries can be directed to the corresponding authors. The computer R codes for the processing and analysis of this study are available at https://github.com/Hao-Zou-lab/Stemness.Subtype.
Tumor mutation load
Immune checkpoint blockade
Cancer stem cells
Single-cell RNA sequencing
Machine learning algorithms
Androgen deprivation therapy
One-class logistic regression
The Cancer Genome Atlas-Prostate Adenocarcinoma
Pathological T stage
Somatic copy number alteration
Fragment genomic alteration
- K–M analysis:
Tumor immune microenvironment
High immune infiltration;
Low immune infiltration
West Coast Dream Team
Copy number alterations
Gene set variation analysis
Gene set enrichment analysis
Ingenuity pathway analysis
Area under the receiver operating characteristic curve
Weighted gene co-expression network analysis
Stemness-related risk scores
Pluripotent stem cells
Gene Expression Omnibus
Progenitor Cell Biology Consortium
Database of Genotypes and Phenotypes
National Center for Biotechnology Information
European Genome-Phenome Archive
Benign prostatic hyperplasia
Cumulative distribution function
Proportion of ambiguous clustering algorithm
Topological overlap matrix
Miller K, Ortiz A, Pinheiro P, Bandi P, Minihan A, Fuchs H, Martinez Tyson D, Tortolero-Luna G, Fedewa S, Jemal A, et al. Cancer statistics for the US Hispanic/Latino population, 2021. Cancer J Clin. 2021;71(6):466–87.
Siegel R, Miller K, Fuchs H, Jemal A. Cancer statistics, 2022. Cancer J Clin. 2022;72(1):7–33.
Bray F, Ferlay J, Soerjomataram I, Siegel R, Torre L, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Cancer J Clin. 2018;68(6):394–424.
Chen W, Zheng R, Baade P, Zhang S, Zeng H, Bray F, Jemal A, Yu X, He J. Cancer statistics in China, 2015. Cancer J Clin. 2016;66(2):115–32.
Marhold M, Kramer G, Krainer M, Le Magnen C. The prostate cancer landscape in Europe: current challenges, future opportunities. Cancer Lett. 2022;526:304–10.
Zheng K, Chen S, Hu X. Peroxisome proliferator-activated receptor gamma coactivator-1 alpha: a double-edged sword in prostate cancer. Curr Cancer Drug Targets. 2022;22(7):541–59.
Topalian S, Hodi F, Brahmer J, Gettinger S, Smith D, McDermott D, Powderly J, Sosman J, Atkins M, Leming P, et al. Five-year survival and correlates among patients with advanced melanoma, renal cell carcinoma, or non-small cell lung cancer treated with nivolumab. JAMA Oncol. 2019;5(10):1411–20.
Motzer RJ, Penkov K, Haanen J, Rini B, Albiges L, Campbell MT, Venugopal B, Kollmannsberger C, Negrier S, Uemura M, et al. Avelumab plus axitinib versus sunitinib for advanced renal-cell carcinoma. N Engl J Med. 2019;380(12):1103–15.
Runcie K, Dallos M. Prostate cancer immunotherapy-finally in from the cold? Curr Oncol Rep. 2021;23(8):88.
Rebuzzi SE, Rescigno P, Catalano F, Mollica V, Vogl UM, Marandino L, Massari F, Pereira Mestre R, Zanardi E, Signori A, et al. Immune checkpoint inhibitors in advanced prostate cancer: current data and future perspectives. Cancers. 2022;14(5):1245.
Zheng K, Gao L, Hao J, Zou X, Hu X. An immunotherapy response prediction model derived from proliferative CD4(+) T cells and antigen-presenting monocytes in ccRCC. Front Immunol. 2022;13:972227.
Qian Y, Feng D, Wang J, Wei W, Wei Q, Han P, Yang L. Establishment of cancer-associated fibroblasts-related subtypes and prognostic index for prostate cancer through single-cell and bulk RNA transcriptome. Sci Rep. 2023;13(1):9016.
Wen XY, Wang RY, Yu B, Yang Y, Yang J, Zhang HC. Integrating single-cell and bulk RNA sequencing to predict prognosis and immunotherapy response in prostate cancer. Sci Rep. 2023;13(1):15597.
Merriel S, Pocock L, Gilbert E, Creavin S, Walter F, Spencer A, Hamilton W. Systematic review and meta-analysis of the diagnostic accuracy of prostate-specific antigen (PSA) for the detection of prostate cancer in symptomatic patients. BMC Med. 2022;20(1):54.
Coleman IM, DeSarkar N, Morrissey C, Xin L, Roudier MP, Sayar E, Li D, Corey E, Haffner MC, Nelson PS. Therapeutic implications for intrinsic phenotype classification of metastatic castration-resistant prostate cancer. Clin Cancer Res. 2022;28(14):3127–40.
Dal Pra A, Ghadjar P, Hayoz S, Liu VYT, Spratt DE, Thompson DJS, Davicioni E, Huang HC, Zhao X, Liu Y, et al. Validation of the Decipher genomic classifier in patients receiving salvage radiotherapy without hormone therapy after radical prostatectomy—an ancillary study of the SAKK 09/10 randomized clinical trial. Ann Oncol. 2022;33(9):950–8.
Cucchiara V, Cooperberg M, Dall’Era M, Lin D, Montorsi F, Schalken J, Evans C. Genomic markers in prostate cancer decision making. Eur Urol. 2018;73(4):572–82.
Zhao SG, Chang SL, Erho N, Yu M, Lehrer J, Alshalalfa M, Speers C, Cooperberg MR, Kim W, Ryan CJ, et al. Associations of luminal and basal subtyping of prostate cancer with prognosis and response to androgen deprivation therapy. JAMA Oncol. 2017;3(12):1663–72.
Gore JL, du Plessis M, Santiago-Jiménez M, Yousefi K, Thompson DJS, Karsh L, Lane BR, Franks M, Chen DYT, Bandyk M, et al. Decipher test impacts decision making among patients considering adjuvant and salvage treatment after radical prostatectomy: interim results from the multicenter prospective PRO-IMPACT study. Cancer. 2017;123(15):2850–9.
Shackleton M, Quintana E, Fearon ER, Morrison SJ. Heterogeneity in cancer: cancer stem cells versus clonal evolution. Cell. 2009;138(5):822–9.
Wolf I, Gratzke C, Wolf P. Prostate cancer stem cells: clinical aspects and targeted therapies. Front Oncol. 2022;12:935715.
Escudero-Lourdes C, Alvarado-Morales I, Tokar EJ. Stem cells as target for prostate cancer therapy: opportunities and challenges. Stem Cell Rev Rep. 2022;18(8):2833–51.
Visvader JE. Cells of origin in cancer. Nature. 2011;469(7330):314–22.
Zhang D, Tang DG, Rycaj K. Cancer stem cells: regulation programs, immunological properties and immunotherapy. Semin Cancer Biol. 2018;52(Pt 2):94–106.
Heidegger I, Fotakis G, Offermann A, Goveia J, Daum S, Salcher S, Noureen A, Timmer-Bosscha H, Schäfer G, Walenkamp A, et al. Comprehensive characterization of the prostate tumor microenvironment identifies CXCR4/CXCL12 crosstalk as a novel antiangiogenic therapeutic target in prostate cancer. Mol Cancer. 2022;21(1):132.
Wong HY, Sheng Q, Hesterberg AB, Croessmann S, Rios BL, Giri K, Jackson J, Miranda AX, Watkins E, Schaffer KR, et al. Single cell analysis of cribriform prostate cancer reveals cell intrinsic and tumor microenvironmental pathways of aggressive disease. Nat Commun. 2022;13(1):6036.
Song H, Weinstein HNW, Allegakoen P, Wadsworth MH 2, Xie J, Yang H, Castro EA, Lu KL, Stohr BA, Feng FY, et al. Single-cell analysis of human primary prostate cancer reveals the heterogeneity of tumor-associated epithelial cell states. Nat Commun. 2022;13(1):141.
Tuong ZK, Loudon KW, Berry B, Richoz N, Jones J, Tan X, Nguyen Q, George A, Hori S, Field S, et al. Resolving the immune landscape of human prostate at a single-cell level in health and cancer. Cell Rep. 2021;37(12):110132.
Chen S, Zhu G, Yang Y, Wang F, Xiao YT, Zhang N, Bian X, Zhu Y, Yu Y, Liu F, et al. Single-cell analysis reveals transcriptomic remodellings in distinct cell types that contribute to human prostate cancer progression. Nat Cell Biol. 2021;23(1):87–98.
Gulati GS, Sikandar SS, Wesche DJ, Manjunath A, Bharadwaj A, Berger MJ, Ilagan F, Kuo AH, Hsieh RW, Cai S, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science. 2020;367(6476):405–11.
Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, Kamińska B, Huelsken J, Omberg L, Gevaert O, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell. 2018;173(2):338–54.e315.
Taylor BS, Schultz N, Hieronymus H, Gopalan A, Xiao Y, Carver BS, Arora VK, Kaushik P, Cerami E, Reva B, et al. Integrative genomic profiling of human prostate cancer. Cancer Cell. 2010;18(1):11–22.
Luca BA, Brewer DS, Edwards DR, Edwards S, Whitaker HC, Merson S, Dennis N, Cooper RA, Hazell S, Warren AY, et al. DESNT: a poor prognosis category of human prostate cancer. Eur Urol Focus. 2018;4(6):842–50.
Grasso CS, Wu YM, Robinson DR, Cao X, Dhanasekaran SM, Khan AP, Quist MJ, Jing X, Lonigro RJ, Brenner JC, et al. The mutational landscape of lethal castration-resistant prostate cancer. Nature. 2012;487(7406):239–43.
Friedrich M, Wiedemann K, Reiche K, Puppel SH, Pfeifer G, Zipfel I, Binder S, Köhl U, Müller GA, Engeland K, et al. The role of lncRNAs TAPIR-1 and -2 as diagnostic markers and potential therapeutic targets in prostate cancer. Cancers. 2020;12(5):1122.
Tomlins SA, Mehra R, Rhodes DR, Cao X, Wang L, Dhanasekaran SM, Kalyana-Sundaram S, Wei JT, Rubin MA, Pienta KJ, et al. Integrative molecular concept modeling of prostate cancer progression. Nat Genet. 2007;39(1):41–51.
Gerhauser C, Favero F, Risch T, Simon R, Feuerbach L, Assenov Y, Heckmann D, Sidiropoulos N, Waszak SM, Hübschmann D, et al. Molecular evolution of early-onset prostate cancer identifies molecular risk markers and clinical trajectories. Cancer Cell. 2018;34(6):996–1011.e1018.
Ross-Adams H, Lamb AD, Dunning MJ, Halim S, Lindberg J, Massie CM, Egevad LA, Russell R, Ramos-Montoya A, Vowler SL, et al. Integration of copy number and transcriptomics provides risk stratification in prostate cancer: a discovery and validation cohort study. EBioMedicine. 2015;2(9):1133–44.
Abida W, Cyrta J, Heller G, Prandi D, Armenia J, Coleman I, Cieslik M, Benelli M, Robinson D, Van Allen EM, et al. Genomic correlates of clinical outcome in advanced prostate cancer. Proc Natl Acad Sci USA. 2019;116(23):11428–36.
Quigley DA, Dang HX, Zhao SG, Lloyd P, Aggarwal R, Alumkal JJ, Foye A, Kothari V, Perry MD, Bailey AM, et al. Genomic hallmarks and structural variation in metastatic prostate cancer. Cell. 2018;174(3):758-769.e759.
Truica MI, Burns MC, Han H, Abdulkadir SA. Turning up the heat on MYC: progress in small-molecule inhibitors. Cancer Res. 2021;81(2):248–53.
Deng L, Qian G, Zhang S, Zheng H, Fan S, Lesinski GB, Owonikoko TK, Ramalingam SS, Sun SY. Inhibition of mTOR complex 1/p70 S6 kinase signaling elevates PD-L1 levels in human cancer cells through enhancing protein stabilization accompanied with enhanced β-TrCP degradation. Oncogene. 2019;38(35):6270–82.
Brunet JP, Tamayo P, Golub TR, Mesirov JP. Metagenes and molecular pattern discovery using matrix factorization. Proc Natl Acad Sci USA. 2004;101(12):4164–9.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.
Long Q, Xu J, Osunkoya AO, Sannigrahi S, Johnson BA, Zhou W, Gillespie T, Park JY, Nam RK, Sugar L, et al. Global transcriptome analysis of formalin-fixed prostate cancer specimens identifies biomarkers of disease recurrence. Cancer Res. 2014;74(12):3228–37.
Hellmann MD, Callahan MK, Awad MM, Calvo E, Ascierto PA, Atmaca A, Rizvi NA, Hirsch FR, Selvaggi G, Szustakowski JD, et al. Tumor mutational burden and efficacy of nivolumab monotherapy and in combination with ipilimumab in small-cell lung cancer. Cancer Cell. 2018;33(5):853–61.e854.
Nymoen DA, Holth A, Hetland Falkenthal TE, Tropé CG, Davidson B. CIAPIN1 and ABCA13 are markers of poor survival in metastatic ovarian serous carcinoma. Mol Cancer. 2015;14:44.
Maxwell KN, Cheng HH, Powers J, Gulati R, Ledet EM, Morrison C, Le A, Hausler R, Stopfer J, Hyman S, et al. Inherited TP53 variants and risk of prostate cancer. Eur Urol. 2022;81(3):243–50.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013;14:7.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Ben-Salem S, Venkadakrishnan VB, Heemers HV. Novel insights in cell cycle dysregulation during prostate cancer progression. Endocr Relat Cancer. 2021;28(6):R141–55.
Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021;22(6):bbab260.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, Liu XS. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 2020;12(1):21.
Hoshida Y, Brunet JP, Tamayo P, Golub TR, Mesirov JP. Subclass mapping: identifying common subtypes in independent disease data sets. PLoS ONE. 2007;2(11):e1195.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.
Ji G, He S, Huang C, Gong Y, Li X, Zhou L. Upregulation of ATP binding cassette subfamily C member 5 facilitates prostate cancer progression and enzalutamide resistance via the CDK1-mediated AR Ser81 phosphorylation pathway. Int J Biol Sci. 2021;17(7):1613–28.
Das R, Sjöström M, Shrestha R, Yogodzinski C, Egusa EA, Chesner LN, Chen WS, Chou J, Dang DK, Swinderman JT, et al. An integrated functional and clinical genomics approach reveals genes driving aggressive metastatic prostate cancer. Nat Commun. 2021;12(1):4601.
Mosler T, Baymaz HI, Gräf JF, Mikicic I, Blattner G, Bartlett E, Ostermaier M, Piccinno R, Yang J, Voigt A, et al. PARP1 proximity proteomics reveals interaction partners at stressed replication forks. Nucleic Acids Res. 2022;50(20):11600–18.
Chen L, Song Y, Hou T, Li X, Cheng L, Li Y, Xing Y. Circ_0004087 interaction with SND1 promotes docetaxel resistance in Prostate cancer by boosting the mitosis error correction mechanism. J Exp Clin Cancer Res. 2022;41(1):194.
Labbé DP, Sweeney CJ, Brown M, Galbo P, Rosario S, Wadosky KM, Ku SY, Sjöström M, Alshalalfa M, Erho N, et al. TOP2A and EZH2 provide early detection of an aggressive prostate cancer subgroup. Clin Cancer Res. 2017;23(22):7072–83.
Tamura K, Furihata M, Tsunoda T, Ashida S, Takata R, Obara W, Yoshioka H, Daigo Y, Nasu Y, Kumon H, et al. Molecular features of hormone-refractory prostate cancer cells by genome-wide gene expression profiles. Cancer Res. 2007;67(11):5117–25.
Sharma NL, Massie CE, Ramos-Montoya A, Zecchini V, Scott HE, Lamb AD, MacArthur S, Stark R, Warren AY, Mills IG, et al. The androgen receptor induces a distinct transcriptional program in castration-resistant prostate cancer in man. Cancer Cell. 2013;23(1):35–47.
Alumkal JJ, Sun D, Lu E, Beer TM, Thomas GV, Latour E, Aggarwal R, Cetnar J, Ryan CJ, Tabatabaei S, et al. Transcriptional profiling identifies an androgen receptor activity-low, stemness program associated with enzalutamide resistance. Proc Natl Acad Sci USA. 2020;117(22):12315–23.
Jain S, Lyons CA, Walker SM, McQuaid S, Hynes SO, Mitchell DM, Pang B, Logan GE, McCavigan AM, O’Rourke D, et al. Validation of a metastatic assay using biopsies to improve risk stratification in patients with prostate cancer treated with radical radiation therapy. Ann Oncol. 2018;29(1):215–22.
Erho N, Crisan A, Vergara IA, Mitra AP, Ghadessi M, Buerki C, Bergstralh EJ, Kollmeyer T, Fink S, Haddad Z, et al. Discovery and validation of a prostate cancer genomic classifier that predicts early metastasis following radical prostatectomy. PLoS ONE. 2013;8(6):e66855.
Mortensen MM, Høyer S, Lynnerup AS, Ørntoft TF, Sørensen KD, Borre M, Dyrskjøt L. Expression profiling of prostate cancer tissue delineates genes associated with recurrence after prostatectomy. Sci Rep. 2015;5:16018.
Fraser M, Sabelnykova VY, Yamaguchi TN, Heisler LE, Livingstone J, Huang V, Shiah YJ, Yousif F, Lin X, Masella AP, et al. Genomic hallmarks of localized, non-indolent prostate cancer. Nature. 2017;541(7637):359–64.
Nabila E. Expression array for Multi-omics molecular profiling of primary prostate adenocarcinoma. BioStudies. E-MTAB-6128. 2018. https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-6128.
Motzer RJ, Robbins PB, Powles T, Albiges L, Haanen JB, Larkin J, Mu XJ, Ching KA, Uemura M, Pal SK, et al. Avelumab plus axitinib versus sunitinib in advanced renal cell carcinoma: biomarker analysis of the phase 3 JAVELIN renal 101 trial. Nat Med. 2020;26(11):1733–41.
Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, Liu XQ, Sher X, Jung H, Lee M, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. 2018;24(9):1449–58.
Snyder A, Nathanson T, Funt SA, Ahuja A, Buros Novik J, Hellmann MD, Chang E, Aksoy BA, Al-Ahmadie H, Yusko E, et al. Contribution of systemic and somatic factors to clinical response and resistance to PD-L1 blockade in urothelial cancer: an exploratory multi-omic analysis. PLoS Med. 2017;14(5):e1002309.
Gide TN, Quek C, Menzies AM, Tasker AT, Shang P, Holst J, Madore J, Lim SY, Velickovic R, Wongchenko M, et al. Distinct Immune cell populations define response to anti-PD-1 monotherapy and anti-PD-1/anti-CTLA-4 combined therapy. Cancer Cell. 2019;35(2):238–55.e236.
Jung H, Kim HS, Kim JY, Sun JM, Ahn JS, Ahn MJ, Park K, Esteller M, Lee SH, Choi JK. DNA methylation loss promotes immune evasion of tumours with high mutation and copy number load. Nat Commun. 2019;10(1):4278.
Pusztai L, Yau C, Wolf DM, Han HS, Du L, Wallace AM, String-Reasor E, Boughey JC, Chien AJ, Elias AD, et al. Durvalumab with olaparib and paclitaxel for high-risk HER2-negative stage II/III breast cancer: results from the adaptively randomized I-SPY2 trial. Cancer Cell. 2021;39(7):989–98.e985.
Cho JW, Hong MH, Ha SJ, Kim YJ, Cho BC, Lee I, Kim HR. Genome-wide identification of differentially methylated promoters and enhancers associated with response to anti-PD-1 therapy in non-small cell lung cancer. Exp Mol Med. 2020;52(9):1550–63.
Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell. 2016;165(1):35–44.
Miao D, Margolis CA, Gao W, Voss MH, Li W, Martini DJ, Norton C, Bossé D, Wankowicz SM, Cullen D, et al. Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Science. 2018;359(6377):801–6.
Zhao J, Chen AX, Gartrell RD, Silverman AM, Aparicio L, Chu T, Bordbar D, Shan D, Samanamud J, Mahajan A, et al. Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma. Nat Med. 2019;25(3):462–9.
Phillips D, Matusiak M, Gutierrez BR, Bhate SS, Barlow GL, Jiang S, Demeter J, Smythe KS, Pierce RH, Fling SP, et al. Immune cell topography predicts response to PD-1 blockade in cutaneous T cell lymphoma. Nat Commun. 2021;12(1):6726.
Choueiri TK, Fishman MN, Escudier B, McDermott DF, Drake CG, Kluger H, Stadler WM, Perez-Gracia JL, McNeel DG, Curti B, et al. Immunomodulatory activity of nivolumab in metastatic renal cell carcinoma. Clin Cancer Res. 2016;22(22):5461–71.
Obradovic A, Graves D, Korrer M, Wang Y, Roy S, Naveed A, Xu Y, Luginbuhl A, Curry J, Gibson M, et al. Immunostimulatory cancer-associated fibroblast subpopulations can predict immunotherapy response in head and neck cancer. Clin Cancer Res. 2022;28(10):2094–109.
Pomeranz Krummel DA, Nasti TH, Izar B, Press RH, Xu M, Lowder L, Kallay L, Rupji M, Rosen H, Su J, et al. Impact of sequencing radiation therapy and immune checkpoint inhibitors in the treatment of melanoma brain metastases. Int J Radiat Oncol Biol Phys. 2020;108(1):157–63.
Chiappinelli KB, Strissel PL, Desrichard A, Li H, Henke C, Akman B, Hein A, Rote NS, Cope LM, Snyder A, et al. Inhibiting DNA methylation causes an interferon response in cancer via dsRNA including endogenous retroviruses. Cell. 2017;169(2):361.
Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-Arnon L, Zimmer L, Gutzmer R, Satzger I, Loquai C, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med. 2019;25(12):1916–27.
Garcia-Diaz A, Shin DS, Moreno BH, Saco J, Escuin-Ordinas H, Rodriguez GA, Zaretsky JM, Sun L, Hugo W, Wang X, et al. Interferon receptor signaling pathways regulating PD-L1 and PD-L2 expression. Cell Rep. 2017;19(6):1189–201.
Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, Ross-Macdonald P, Berger AC, Jegede OA, Elagina L, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med. 2020;26(6):909–18.
Ascierto ML, McMiller TL, Berger AE, Danilova L, Anders RA, Netto GJ, Xu H, Pritchard TS, Fan J, Cheadle C, et al. The intratumoral balance between metabolic and immunologic gene expression is associated with Anti-PD-1 response in patients with renal cell carcinoma. Cancer Immunol Res. 2016;4(9):726–33.
Patil NS, Nabet BY, Müller S, Koeppen H, Zou W, Giltnane J, Au-Yeung A, Srivats S, Cheng JH, Takahashi C, et al. Intratumoral plasma cells predict outcomes to PD-L1 blockade in non-small cell lung cancer. Cancer Cell. 2022;40(3):289-300.e284.
Keenan TE, Guerriero JL, Barroso-Sousa R, Li T, O’Meara T, Giobbie-Hurder A, Tayob N, Hu J, Severgnini M, Agudo J, et al. Molecular correlates of response to eribulin and pembrolizumab in hormone receptor-positive metastatic breast cancer. Nat Commun. 2021;12(1):5563.
Motzer RJ, Banchereau R, Hamidi H, Powles T, McDermott D, Atkins MB, Escudier B, Liu LF, Leng N, Abbas AR, et al. Molecular subsets in renal cancer determine outcome to checkpoint and angiogenesis blockade. Cancer Cell. 2020;38(6):803-817.e804.
Lauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, Salim M, Vallon-Christersson J, Törngren T, Kvist A, et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat Commun. 2017;8(1):1738.
van den Ende T, de Clercq NC, van Berge Henegouwen MI, Gisbertz SS, Geijsen ED, Verhoeven RHA, Meijer SL, Schokker S, Dings MPG, Bergman J, et al. Neoadjuvant chemoradiotherapy combined with atezolizumab for resectable esophageal adenocarcinoma: a single-arm phase II feasibility trial (PERFECT). Clin Cancer Res. 2021;27(12):3351–9.
Motzer RJ, Rini BI, McDermott DF, Redman BG, Kuzel TM, Harrison MR, Vaishampayan UN, Drabkin HA, George S, Logan TF, et al. Nivolumab for metastatic renal cell carcinoma: results of a randomized phase II trial. J Clin Oncol. 2015;33(13):1430–7.
Motzer RJ, Tannir NM, McDermott DF, Arén Frontera O, Melichar B, Choueiri TK, Plimack ER, Barthélémy P, Porta C, George S, et al. Nivolumab plus ipilimumab versus sunitinib in advanced renal-cell carcinoma. N Engl J Med. 2018;378(14):1277–90.
Kraehenbuehl L, Holland A, Armstrong E, O’Shea S, Mangarin L, Chekalil S, Johnston A, Bomalaski JS, Erinjeri JP, Barker CA, et al. Pilot trial of arginine deprivation plus nivolumab and ipilimumab in patients with metastatic uveal melanoma. Cancers. 2022;14(11):2638.
Wolf DM, Yau C, Wulfkuhle J, Brown-Swigart L, Gallagher RI, Lee PRE, Zhu Z, Magbanua MJ, Sayaman R, O’Grady N, et al. Redefining breast cancer subtypes to guide treatment prioritization and maximize response: predictive biomarkers across 10 cancer therapies. Cancer Cell. 2022;40(6):609–23.e606.
Liu S, Knochelmann HM, Lomeli SH, Hong A, Richardson M, Yang Z, Lim RJ, Wang Y, Dumitras C, Krysan K, et al. Response and recurrence correlates in individuals treated with neoadjuvant anti-PD-1 therapy for resectable oral cavity squamous cell carcinoma. Cell Rep Med. 2021;2(10):100411.
Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, Tian T, Wei Z, Madan S, Sullivan RJ, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat Med. 2018;24(10):1545–9.
Mahoney KM, Ross-Macdonald P, Yuan L, Song L, Veras E, Wind-Rotolo M, McDermott DF, Stephen Hodi F, Choueiri TK, Freeman GJ. Soluble PD-L1 as an early marker of progressive disease on nivolumab. J Immunother Cancer. 2022;10(2):e003527.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8.
Riaz N, Havel JJ, Makarov V, Desrichard A, Urba WJ, Sims JS, Hodi FS, Martín-Algarra S, Mandal R, Sharfman WH, et al. Tumor and microenvironment evolution during immunotherapy with nivolumab. Cell. 2017;171(4):934-949.e916.
Pan KF, Lee WJ, Chou CC, Yang YC, Chang YC, Chien MH, Hsiao M, Hua KT. Direct interaction of β-catenin with nuclear ESM1 supports stemness of metastatic prostate cancer. Embo J. 2021;40(4):e105450.
Shorning BY, Dass MS, Smalley MJ, Pearson HB. The PI3K-AKT-mTOR pathway and prostate cancer: at the crossroads of AR, MAPK, and WNT signaling. Int J Mol Sci. 2020;21(12):4507.
Horning AM, Wang Y, Lin CK, Louie AD, Jadhav RR, Hung CN, Wang CM, Lin CL, Kirma NB, Liss MA, et al. Single-cell RNA-seq reveals a subpopulation of prostate cancer cells with enhanced cell-cycle-related transcription and attenuated androgen response. Cancer Res. 2018;78(4):853–64.
Chen K, Zhang C, Ling S, Wei R, Wang J, Xu X. The metabolic flexibility of quiescent CSC: implications for chemotherapy resistance. Cell Death Dis. 2021;12(9):835.
Matthews HK, Bertoli C, de Bruin RAM. Cell cycle control in cancer. Nat Rev Mol Cell Biol. 2022;23(1):74–88.
Gupta N, Huang TT, Horibata S, Lee JM. Cell cycle checkpoints and beyond: exploiting the ATR/CHK1/WEE1 pathway for the treatment of PARP inhibitor-resistant cancer. Pharmacol Res. 2022;178:106162.
de Morree A, Rando TA. Regulation of adult stem cell quiescence and its functions in the maintenance of tissue integrity. Nat Rev Mol Cell Biol. 2023;24(5):334–54.
Tsao T, Beretov J, Ni J, Bai X, Bucci J, Graham P, Li Y. Cancer stem cells in prostate cancer radioresistance. Cancer Lett. 2019;465:94–104.
Dunning NL, Laversin SA, Miles AK, Rees RC. Immunotherapy of prostate cancer: should we be targeting stem cells and EMT? Cancer Immunol Immunother. 2011;60(8):1181–93.
Cha HR, Lee JH, Ponnazhagan S. Revisiting immunotherapy: a focus on prostate cancer. Cancer Res. 2020;80(8):1615–23.
Chen F, Wendl MC, Wyczalkowski MA, Bailey MH, Li Y, Ding L. Moving pan-cancer studies from basic research toward the clinic. Nat Cancer. 2021;2(9):879–90.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005;102(43):15545–50.
Szklarczyk D, Kirsch R, Koutrouli M, Nastou K, Mehryary F, Hachilif R, Gable AL, Fang T, Doncheva NT, Pyysalo S, et al. The STRING database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 2023;51(D1):D638–46.
Șenbabaoğlu Y, Michailidis G, Li JZ. Critical limitations of consensus clustering in class discovery. Sci Rep. 2014;4:6207.
The authors are grateful to the Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), UCSC Xena, cBioPortal, ArrayExpress, European Genome-phenome Archive (EGA), National Center for Biotechnology Information (NCBI), Database of Genotypes and Phenotypes (dbGaP), Progenitor Cell Biology Consortium (PCBC), and pioneer researchers who provided the datasets. The part computations in this study were run on the Siyuan-1 cluster, supported by the Center for High Performance Computing at Shanghai Jiao Tong University.
This work was supported in part by the National Natural Science Foundation of China [82170045 to JH], Innovative Research Team of High-level Local Universities in Shanghai [SHSMU-ZLCX20212301 to JH], and Shanghai Municipal Health Commission Project [202040375 to XH].
Ethics approval and consent to participate
The institutional review board of Shanghai Sixth People’s Hospital approved this study, and informed consent was obtained from all patients.
Consent for publication
Written informed consent for publication of their clinical details and/or clinical images was obtained from the patient/parent/guardian/relative of the patient.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Characteristics of scRNA-seq datasets and bulk RNA-seq datasets enrolled in this study. Table S2. Details of 9-gene stemness subtype predictor. Figure S1. The holistic design of the current study. Figure S2. Uniform manifold approximation and projection (UMAP) of all cells captured across 5 single-cell RNA-seq (scRNA-seq) datasets. Figure S3. Correlation of stemness levels with clinical, pathological, and molecular features in patients with prostate cancer (PCA). Figure S4. Correlation of stemness levels with tumor immune microenvironment (TIME) in PCA patients. Figure S5. Development of 18 signatures for stemness classification. Figure S6. Identification of three PCA stemness subtypes based on stemness signatures. Figure S7. Comparison of molecular features among three PCA stemness subtypes. Figure S8. Signaling pathways and functional annotations of three PCA stemness subtypes. Figure S9. Comparison of drug sensitivities and TIME patterns among three PCA stem subtypes. Figure S10. Construction of stemness subtype predictor. Figure S11. Validation of stemness subtype predictor. Figure S12. Conservation of the three stemness subtypes in pan-tumors. Figure S13. The impact of three pan-cancer stemness subtypes on ICB treatment and their predictive performance for ICB efficacy. Additional Remarks A-G. Methods S1.
Eighteen stemness signatures for stemness classification. Data S2. Prediction performances of 9- gene stemness subtype predictor based on 100 machine learning algorithms (MLs). Data S3. Details of 100 machine learning algorithms.
About this article
Cite this article
Zheng, K., Hai, Y., Xi, Y. et al. Integrative multi-omics analysis unveils stemness-associated molecular subtypes in prostate cancer and pan-cancer: prognostic and therapeutic significance. J Transl Med 21, 789 (2023). https://doi.org/10.1186/s12967-023-04683-6