The RNA-seq expression profiles, mutation data, and corresponding clinical information of PCa were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database, which contained 499 PCa specimens and 52 normal tissue specimens. Characteristics of the PCa patients from the TCGA data were listed in Table 1. Furthermore, from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database, we downloaded gene expression and complete biochemical recurrence (BCR) data for 248 PCa samples in GSE116918 (Jain et al., 2018). Then we validated our results in an independent external validation dataset that combined the transcriptional expression of key prognostic genes and clinical data. We acquired a total of 1,011 protein-coding EMT-related genes from dbEMT 2.0 (https://dbemt.bioinfo-minzhao.org/index.html) (Zhao et al., 2019).

Differentially expressed genes related to EMT from TCGA datasets were identified between PCa and normal prostate samples using the ‘limma’ R package. The genes were considered significant if they had an adjusted P-value <0.05 and log |fold change (FC)| > 1. The functional enrichment analysis of the differentially expressed EMT-related genes was systematically explored using Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses with the ‘clusterProfiler’ R package (Yu et al., 2012).

The search tool for the retrieval of interacting genes/proteins (STRING) database (http://string-db.org) was used to construct protein-protein interactions (PPI) among the differentially expressed EMT-related genes [18]. The most important module of the PPI network was then identified using the Molecular Complex Detection (MCODE) plugin of Cytoscape with default parameters [19].

Non-negative matrix factorization (NMF) was used to classify patients into different EMT regulator patterns based on the gene expression profiling of differentially expressed EMT-related genes. The ‘NMF’ R package was utilized for this purpose (Gaujoux and Seoighe, 2010). The number of clusters (k) ranged from 2 to 10, and the optimal number (k = 2) was determined based on factorization ranking parameters such as cophenetic and dispersions as well as consensus heatmap. The study analyzed the differences in progression-free survival (PFS) between distinct clusters. Additionally, the immune infiltration landscape of different clusters was elucidated using the single-sample gene set enrichment analysis (ssGSEA) algorithm (Subramanian et al., 2005). The gene expression and the clinicopathological distribution of the clusters were also visualized using the heatmaps. Finally, the expression level of PD-1 was evaluated across different PCa patterns.

Differentially expressed genes (DEGs) were identified between the different EMT regulator patterns, with an adjusted P-value <0.05 and log |fold change (FC)| > 2. The genes related to the EMT regulator pattern were then analyzed using univariate Cox regression analysis to determine their association with PFS. Additionally, patients were divided into different gene clusters based on the expression of prognostic DEGs using the NMF unsupervised clustering approach.

The EMT prognostic signature (EPS) was constructed using the TCGA dataset. To minimize over-fitting prognostic characteristics, the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was employed through the ‘glmnet’ R package. The study identified key prognostic genes through multivariate Cox analysis and constructed an EPS using the formula: risk score = ∑ (Coefi*Expi), where Coefi and Expi represent the regression coefficient and corresponding expression value of each gene. According to the median risk score, 595 patients in the TCGA dataset were sorted into low-risk group and high-risk group, and then subjected to survival analysis. Furthermore, the performance of the constructed EPS was tested using time-dependent receiver operating characteristic (ROC) curve analysis, principal component analysis (PCA), and t-distributed stochastic neighbor embedding (t-SNE) with the ‘timeROC’ and ‘Rtsne’ R packages. The EPS’s utility was validated using the GSE116918 dataset as an independent external cohort. With the 'survminer’ R program, an optimal cut-off value for survival analysis was established.

To determine the practical utility for the EPS, correlation analyses, and stratified survival analyses in the TCGA dataset were carried out. Next, to ascertain if the signature may function as an independent prognostic factor, we conducted univariate and multivariate Cox regression analyses. A clinical nomogram was developed using the 'rms’ R package to predict PFS by combining independent prognostic factors. The precision of the nomogram prediction outcomes was evaluated with calibration curves as well as decision curve analysis (DCA) via the 'ggDCA’ R package.

The tumor mutation burden (TMB) value of each PCa sample was calculated with the 'maftools’ R package (Mayakonda et al., 2018) after downloading somatic variant data. Additionally, we conducted differential and survival analyses of TMB between two risk groups. We also utilized Pearson’s method to examine the correlation between the risk score and various genomic heterogeneities, including homologous recombination defect (HRD), loss of heterozygosity (LOH), microsatellite instability (MSI), purity, mutant-allele tumor heterogeneity (MATH), neoantigens, and ploidy (Bonneville et al., 2017; Thorsson et al., 2018).

Immune infiltration scores were assessed by currently accepted methods, such as TIMER (Li et al., 2020), CIBERSORT (Chen et al., 2018), CIBERSORT-ABS (Tamminga et al., 2020), QUANTISEQ (Chen et al., 2018), MCPCOUNT (Dienstmann et al., 2019), XCELL (Aran et al., 2017), and EPIC (Racle et al., 2017). The relative abundance of immune cells between different risk groups of PCa patients was compared by a ssGSEA. Additionally, the relative contents of 22 tumor-infiltrating immune cells (TIICs) were obtained with the CIBERSORT algorithm. We then conducted correlation analyses between the risk scores and the TIIC contents. Finally, we analyzed the relationship between the two risk groups and the immune subtypes reported by Thorsson et al. This study conducted an immunogenomic analysis of more than a thousand tumor samples covering 33 different cancer types and defined six immunological categories (C1-C6) (Thorsson et al., 2018).

We compared the expression levels of 46 immune checkpoints between two risk groups. Then, we utilized the Cancer Immunome Atlas (TCIA, https://tcia.at/home) database (Charoentong et al., 2017) to quantitatively measure tumor immunogenicity using the immunophenoscore (IPS), which ranges from 0–10. To identify the relationship between drug sensitivity and risk score, we used the “oncoPredict” R package to predict the five most common drug sensitivities (Maeser et al., 2021). The difference in sensitivity score was then compared between the two groups with Wilcoxon signed-rank tests. Subsequently, the anti-cancer drug targets were extracted from the DrugBank (https://go.drugbank.com/) database (Wishart et al., 2018) and expression differences were also analyzed between the two groups. The association of the key prognostic gene with chemotherapeutic drug sensitivity was explored to validate the efficacy of these genes in predicting drug sensitivity using the Gene Set Cancer Analysis (GSCA, https://guolab.wchscu.cn/GSCA) database (Liu et al., 2023).

The protein expression levels of the key prognostic genes in PCa relative to normal tissue were measured via immunohistochemistry (IHC) staining in the Human Protein Atlas database (HPA, https://www.proteinatlas.org/) (Asplund et al., 2012). The IHC images were downloaded from this database. The IHC images were analyzed by calculating the percentage of the reaction areas using ImageJ software, version 1.53k, (Wayne Rasband, the National Institutes of Health in the USA). Data visualization was conducted with the GraphPad Prism software, version 9.0.0.

Additionally, tissues from ten PCa patients were obtained for quantitative real-time polymerase chain reaction (qRT-PCR) assays. This study was approved by the Ethics Committee of The Sixth Affiliated Hospital of Sun Yat-sen University (the Ethical Approval Number: 2022ZSLYEC-468). All patients signed written informed consent forms to use their histopathological samples for research purposes, and the study was conducted in accordance with the Declaration of Helsinki. Total RNA was extracted from cells using TRIzol™ Reagent (Invitrogen, USA) according to the manufacturer’s instructions. cDNA was synthesized from 1 μg of total RNA using a reverse transcription kit (please specify the kit and manufacturer, if available). Quantitative real-time PCR (qRT-PCR) was performed using SYBR GreenER™ qPCR SuperMix (Invitrogen, USA) on a Bio-Rad iCycler system (Bio-Rad Laboratories, CA, USA). The thermal cycling conditions were as follows: initial denaturation at 95°C for 1 min, followed by 35 cycles of denaturation at 95°C for 90 s, annealing at 60°C for 30 s, and extension at 72°C for 30 s. A final extension was performed at 72°C for 10 min. The relative mRNA levels were standardized using the 2^−ΔΔCT technique and compared to GAPDH. The qRT-PCR primer sequences can be found in Table 2.

Information on EMT-related genes in homo sapiens was obtained from the dbEMT2.0 database, with 1,184 genes identified. Of these, 85.4% (n = 1,011) were retained for subsequent analysis. Transcriptome data for these genes were obtained from the TCGA database. Differential expression analysis identified 186 EMT-related DEGs, with 42 upregulated DEGs and 144 downregulated DEGs in PCa relative to normal tissue (Figures 1A,B).

Functional enrichment analysis was utilized to investigate the molecular mechanisms related to the EMT-related DEGs. The biological processes were associated with gland development, cell-substrate adhesion, and mesenchyme development. The cellular component was found to be involved in the collagen-containing extracellular matrix, focal adhesion, and cell-substrate junction. The molecular function was related to receptor-ligand activity, cytokine receptor binding, and transmembrane receptor protein kinase activity (Figure 1C). Leading edge analysis of 186 differentially expressed EMT-related genes was shown in Figure 1D. Additionally, KEGG enrichment analysis revealed the enrichment of several malignancy-associated signaling pathways (Figure 1E).

The PPI networks of 186 EMT-related DEGs were reconstructed using the STRING database. The resulting network of interaction relationships consisted of 161 nodes and 917 edges (with a minimum required interaction score >0.4), indicating the complexity of EMT in PCa progression (Figure 1F). The most important module was identified using the Cytotype parameter Settings with the MCODE plug-in (Figure 1G). According to the KEGG enrichment analysis, the genes in this module were found to be highly correlated with PCa (Figure 1H). This implies that these genes may potentially contribute significantly to the progression of PCa.

Consensus clustering was conducted using the NMF algorithm to differentiate PCa patients with different EMT regulation patterns based on the expression of 186 EMT-related DEGs. The NMF rank survey, including cophenetic and dispersion (Figure 2A), and the consensus matrix heatmap (k = 2–5, Figure 2B), showed that k = 2 was the optimal parameter to categorize the PCa patients into two distinct regulation patterns. There were 309 cases in EMT cluster 1 and 186 cases in EMT cluster 2. Characteristics of cluster 1 and cluster 2 from TCGA data were listed in Table 1. Prognostic analysis revealed that EMT cluster 1 had a more significant survival advantage than EMT cluster 2 (Figure 2C). The heatmap of the expression of EMT-related DEGs demonstrated the distinction between the two clusters (Figure 2D). Next, we utilized the ssGSEA method to analyze the differences in immune cell composition between the two EMT clusters. The study results indicate that the abundance of immune cells, including mast cells, natural killer T cells, neutrophils, regulatory T cells, and follicular/Type1 helper T cells, was significantly higher in EMT cluster 1. On the other hand, the abundance of Gamma-delta (γδ) T cells and Type2 helper T cells was significantly higher in EMT cluster 2 (Figure 2E). Additionally, tumors in EMT cluster 1 showed increased expression of PD-1 and PD-L1 (Figure 2F). These results indicate that two EMT regulatory patterns are distinct in terms of clinical features and immune infiltration levels, suggesting that PCa cases can be classified according to EMT-related genes.

To demonstrate the heterogeneity of different EMT regulator patterns, we compared gene expression profiles between the two clusters and identified 2,512 genes associated with the EMT regulator pattern (Supplementary Figure S1A). Subsequently, a univariate Cox regression was conducted on these genes and identified 873 genes significantly associated with PFS (Supplementary Figure S1B). Consistent with the clustering of EMT regulation patterns, NMF cluster analysis based on the expression of these 873 genes explicitly categorized PCa patients into four clusters, which we termed EMT gene clusters 1-4, respectively (Supplementary Figures 1C, D ). Of the 495 patients with PCa, 84 were clustered into gene cluster 1, which was linked to the best prognosis among the four gene clusters, while patients in gene cluster 2 (n = 178) experienced the worst outcome (Supplementary Figure S1E). The expression profiles of the 873 genes that regulate EMT patterns and their clinical characteristics were illustrated in a heatmap (Supplementary Figure S1F). Significant differences in the expression of EMT-related DEGs were observed in the most important module of the four EMT gene clusters, which is in accordance with the predicted findings of the EMT patterns (Supplementary Figure S1G).

To quantify individual patients’ EMT scores, we constructed an EPS considering the heterogeneity and complexity of PCa patients. We included 873 EMT regulator pattern-related prognostic genes in the LASSO-Cox regression model and screened out nine genes based on the minimum value of λ (Figures 3A,B). A multivariate Cox regression analysis was conducted on nine genes using the Akaike information criterion value. Four genes (MEN1, H2AFZ, UCKL1, and FUS) were selected to establish an EPS (Figure 3C). The risk score was derived with the expression levels and coefficients of the above four genes: risk score = (0.0575 * MEN1) + (0.0150 *H2AFZ) + (0.0765 * UCKL1) + (0.0347 * FUS).

Next, the patients were classified into low-risk and high-risk groups according to their median risk scores. According to the Kaplan-Meier curve, PCa patients in the low-risk group had better PFS (P < 0.001) (Figure 3D). The area under the ROC curve (AUC) values for 1-, 3-, and 5-year survival were 0.786, 0.769, and 0.726, respectively (Figure 3E). The gene expression profiling heatmap showed a significant increase in the expression of four genes in the high-risk stratification (Figure 3F). The plot of survival status and risk score revealed that the high-risk group had a worse PFS rate (Figure 3G). Additionally, we demonstrated the bidirectional distribution of PCa patients in different groups by PCA and t-SNE analyses (Figure 3H).

To verify the prognostic potential of the EPS, a similar analysis was performed on the GSE116918 dataset. Patients were categorized into two groups based on the formula used in the TCGA cohort. The high-risk group exhibited lower survival rates (P < 0.001) (Supplementary Figure S2A). The AUC values of 0.992, 0.873, and 0.613 predicted 1-, 3-, and 5-year BCR, respectively (Supplementary Figure S2B). Additionally, our results were consistent with those of the TCGA dataset in terms of the distribution plots of the expression of the 4 genes, risk score and survival status, PCA, and t-SNE analyses (Supplementary Figures 2C–E). The results indicate that the EPS is a significant biomarker for predicting the survival of PFS and BCR, with potential clinical applications.

To investigate the correlation between the risk scoring system and clinicopathological characteristics, we analyzed differences in various stratified features. The findings revealed that in the TCGA cohort, patients with older age (P < 0.05), higher T-stage (P < 0.01), and higher N-stage (P < 0.001) had higher risk scores (Supplementary Figures 3A, B). Furthermore, patients in gene cluster 2 and cluster 1 had the highest and lowest risk score, which was consistent with their survival outcomes. Next, we investigated the predictive value of risk scores for various clinical characteristics. Furthermore, stratified prognostic analysis revealed that the high-risk group had a worse prognosis in all subgroups except for the N1 stage (Supplementary Figure S3C).

Univariate and multivariate Cox regression analyses were conducted to determine the predictiveness of risk scores. The study found that the T-stage and risk score are independent prognostic indicators of PFS in PCa patients (Figures 4A,B). A nomogram using these independent indicators (T-stage and risk score) was developed to predict PFS (Figure 4C). The calibration curves showed good agreement between the predicted and actual observations of 1-, 3-, and 5-year PFS (Figure 4D). Additionally, the DCA demonstrated that the nomogram provided a greater net benefit compared to other independent factors, indicating its reliability (Figure 4E). To summarize, the predictive nomogram can accurately predict an individual’s survival risk and aid in clinical management.

Since genetic mutations and genomic heterogeneity are the conditions for oncogenesis, we visualized somatic mutations of PCa samples based on the EPS, which showed that SPOP (11%), TP53 (11%), and TTN (10%) had the highest mutation frequencies (Supplementary Figure S4A). Analysis of variance showed that TMB was significantly higher in the high-risk group than in the low-risk group (Supplementary Figure S4B). Spearman’s correlation analysis revealed a positive correlation between TMB and the risk score (Supplementary Figure S4C). The Kaplan-Meier curve also demonstrated that the low-TMB group had a better prognosis (P = 0.01) (Supplementary Figure S4D). By combining TMB and the risk score, PCa patients were categorized into four subgroups for survival assessment. The survival analysis showed that patients with low TMB and low risk had the best prognosis (P < 0.01) (Supplementary Figure S4E), which validates the EPS. Correlation analysis subsequently revealed that the risk score was positively correlated with HRD (R = 0.37, P = 3.5e-16), LOH (R = 0.34, P = 3.3e-14), MSI (R = 0.28, P = 1.7e-10), and purity (R = 0.42, P = 1.3e-07), but not with MATH, neoantigen, and polity (Supplementary Figure S4F). These results may reflect the practical application of the EPS in exploring PCa heterogeneity.

Various levels of immune cell infiltration play a crucial role in the development and progression of cancer. We conducted an analysis of immune cell infiltration using seven algorithms to investigate the correlation between the risk score and the tumor immune microenvironment in patients with PCa (Figure 5A). Besides, the study assessed the abundance of tumor-infiltrating immune cells and immune functions with the ssGSEA algorithm (Figure 5B). Moreover, the results suggest that there may be a difference in immune cell infiltration that leads to altered immune function. The results showed that the low-risk group had significantly higher immune cell and function scores than the high-risk group (Figures 5C,D). Although the stroma score was higher in the low-risk group, there was no significant difference in immune scores between the two subgroups (Figures 5E,F). To further clarify the effect of the risk score on 22 tumor-infiltrating immune cells (TIICs), we calculated the TIIC levels for each sample using the CIBERSORT method. The risk score was shown to be positively associated with Tregs and CD8 T cells, but negatively associated with resting memory CD4 T cells, resting mast cells, and neutrophils (Figure 5G). Upon comparing the relationship between the risk score and the defined immune categories, as per the study conducted by Thorsson et al. (Thorsson et al., 2018), we observed a distribution difference in the proportions of C3 (inflammatory) and C4 (lymphocyte depleted) between the low-risk and high-risk groups (P = 0.001) (Figure 5H). Specifically, C3 was prevalent in the low-risk group, which was significantly associated with good survival, consistent with our findings. C4 was highly clustered in the high-risk group. This result not only reveals the unique features of the PCa immune microenvironment but also adds to previous studies.

Considering the significance of tumor immunotherapy based on immune checkpoint inhibitors (ICIs), we further investigated the differential expression of immune checkpoints between the two groups. In the high-risk group, 13 out of 46 immune checkpoints, including PD-1 (PDCD1) and CTLA4, were upregulated, while 10 immune checkpoints, including PD-L1 (CD274) and PD-L2 (PDCD1LG2), were downregulated (Figure 6A). Additionally, we investigated the immune-related progression-free survival of patients undergoing various treatments. Our findings indicate that in the CTLA4-negative/PD-1-positive subgroups, patients with a low-risk score had a higher IPS (P = 0.03). This higher IPS predicted a more favorable therapeutic outcome (Figure 6B).

Furthermore, to determine the feasibility of using the risk score to personalize chemotherapy for PCa, we explored whether risk scores correlated with IC50 values for five commonly used medications. As shown in Figure 6C, Cisplatin, Docetaxel, Gemcitabine, Paclitaxel, and Vinblastine exhibited greater sensitivity in the low-risk group. From the gene expression data obtained from the DrugBank database, it was observed that nine genes targeted by these drugs had different expression patterns between the two groups (Figure 6D). We next conducted the correlation analysis between risk prognostic-related four key genes and the clinical efficacy of PCa treatments. By analyzing the drug response data from the GDSC database, we noticed that these four genes were significantly linked with sensitivity or resistance to multiple therapeutic drugs and molecular inhibitors using Spearman’s correlation coefficient. The results indicated that the high expression of MEN1, H2AFZ, UCKL1, and FUS produced resistance to 17-AAG (a selective HSP90 inhibitor), and MEN1, H2AFZ, and FUS also produced resistance to TGX221 (a selective PI3K inhibitor) and Trametinib (a selective MEK inhibitor). Simultaneously, MEN1, H2AFZ UCKL1, and FUS were all sensitive to Navitoclax (ABT-263), an effective Bcl-2 inhibitor. MEN1, H2AFZ, and FUS were also sensitive to Tubastatin (a selective HDAC6 inhibitor), Vorinostat (an HDAC inhibitor), and other molecular inhibitors as well as drugs (Figure 6E).

As previously mentioned, the EPS was based on four key prognostic genes, including MEN1, H2AFZ, UCKL1, and FUS. It is noteworthy that MEN1 and H2AFZ are EMT-related genes. We associated the expression levels of these four key genes with clinicopathological features and found that their expression was inversely correlated with T- and N-stage (Figures 7A,B). IHC were introduced to examine the protein expression levels of the four genes in the HPA database, the results revealed a significant increase in their protein expression in PCa tissues compared to normal prostate tissues (all P < 0.05) (Figure 7C). Furthermore, we conducted a real-world experiment to detect the relative mRNA expression levels using qRT-PCR. The study confirmed that the expression levels of MEN1, H2AFZ, UCKL1, and FUS were higher in PCa than in adjacent normal prostate tissues (P < 0.05) (Figure 7D). This finding is consistent with the bioinformatic analyses.