In total, 1,145 patients with prostate adenocarcinoma (PRAD) from TCGA, ICGC, and GEO databases were enrolled in this study. Of these, the data of 540 patients with PRAD from TCGA database were downloaded as the training cohort (https://portal.gdc.cancer.gov/). Two test cohorts were used, namely, test cohort I consisting of 357 patients with PRAD from the ICGC database (https://dcc.icgc.org/projects/PRAD-CA) and test cohort II consisting of 248 patients with PRAD from the GSE116918 database (radical radiotherapy with ADT) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116918). All downloaded data included fragments per kilobase of sequences per million mapped reads (FPKM)-normalized RNA sequencing (RNA-seq) data, clinical characteristics annotation, and follow-up information such as biochemical recurrence (BCR), metastasis (Met), and overall survival (OS). Before further analysis, all RNA-seq data included were log2-transformed and normalized by the R package “sva.”

First, based on the HALLMARKS gene set from the Molecular Signatures Database (MSigDB) (http://www.gsea-msigdb.org/gsea/msigdb/search.jsp), ssGSEA was calculated using the R package “GSVA” in the training set (Lee et al., 2008). The R package “WGCNA” was used to perform WGCNA by using TCGA database mRNA matrix (Langfelder and Horvath, 2008). An adjacency matrix was constructed to describe the correlation strength between the nodes and then transformed to a topological overlap matrix (TOM). Next, hierarchical clustering was performed to identify the modules by setting the minimum number of genes in each module to 60. After merging similar modules, the module genes with the highest correlation to E2F target scores in ssGSEA results were identified. These genes were then combined with the E2F target gene set in HALLMARKS to obtain a new gene set, denoted as gene set A.

Simultaneously, Kaplan–Meier survival and log-rank tests (R package “survival”) were performed using gene set A and BCR information to screen genes associated with prognosis—which were labeled as gene set B. DEGs between TCGA-PRAD tumors and normal tissues were obtained by using R package “limma” (false discovery rate [FDR] < 0.05, log|fold change [FC]|> 1) (Ritchie et al., 2015). E2F target-related candidate prognostic genes were obtained by intersecting DEGs with gene set B. Furthermore, these genes were enrolled for inclusion in LASSO Cox regression (R package “glmnet”) (Tibshirani, 1997). Finally, the remaining genes after 10,000 contractions were selected, and their coefficients were recorded to obtain an E2F target-related gene signature.

All analyses and graphs were executed using R v4.0.3 (http://www.r-project.org) and GraphPad v8.0.3 (https://www.graphpad.com/). The Pearson correlation test was used to calculate the established inter-gene correlations in the gene signature. After the E2F target-related risk score was calculated based on the signature, the training set (TCGA cohort) and the two test sets (cohorts from ICGC and GEO) were divided into high- and low-risk groups according to the median risk score in the training set. Differences between the groups were calculated using t-tests and chi-square tests. Scatter maps and heat maps were used to visualize the survival distribution and gene expression patterns of patients in different groups (R packages “ggplot2” and “scales”). Furthermore, the dimension of high-latitude data was reduced by using t-distributed stochastic neighbor embedding (tSNE) and principal components analysis (PCA) of the R package “Rtsne” to test whether the gene signature could divide patients into different groups (Reich et al., 2008; van der Maaten and Hinton, 2008). The prognostic ability of the gene signature was reflected by the Kaplan–Meier survival analysis, receiver operating characteristic (ROC) curve, and time-dependent ROC (T-ROC) curve (R packages “survival” and “timeROC”). The multivariate Cox stepwise regression models of the R package “survival” were used to test whether the gene signature was an independent predictor of prognosis.

The function between groups was analyzed by GSEA v4.1.0. For this algorithmic analysis, it was verified whether the E2F target pathway was activated in different groups. In addition, GSEA based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases was used to study functional differences according to the DEGs in different groups (R package “GOplot” and “limma”) (Subramanian et al., 2005).

The immune cells were divided into different types based on previous studies. Meanwhile, the differences of immune infiltration were explored by ssGSEA (immune gene set from MSigDB) and the CIBERSORT algorithm (https://cibersortx.stanford.edu/), as well as the immunophenoscore (IPS) and immunophenotyping (Newman et al., 2015; Charoentong et al., 2017; Thorsson et al., 2018). Finally, gene correlation analysis based on the CellMiner database (https://discover.nci.nih.gov/cellminer/home.do) explored the sensitivity between the genes and drugs to screen for new medicines and therapeutic targets (Pommier et al., 2012).

The overall process is shown by a schematic diagram in Figure 1A. In the first step, ssGSEA scores of 50 pathways in each sample were included in univariate Cox regression analysis. For the BCR of the training cohort, the E2F target pathway had the greatest statistical significance (p < 0.0001, Figure 1B). The E2F target pathway score was z-score-transformed and divided into high- and low-z-score groups according to their median value. The incidence of BCR in the high-z-score group was significantly higher than that of the low-group (p < 0.0001, Figure 1C). Moreover, the Kaplan–Meier survival curves revealed a worse prognosis in the high-z-score group (p < 0.001, Figure 1D). The details of the training cohort are presented in Supplementary Table S1.

A total of 35 gene modules were identified by WGCNA, among which the expression of the black module was most related to the increase in the E2F target z-score, with a total of 806 genes (Figures 2A,B and Supplementary Table S2). Then, 1,006 genes related to the E2F target pathway were obtained by combining the 806 genes with those of the E2F target pathway from MSigDB (Supplementary Table S2). After univariate Cox regression analysis, 359 E2F-related prognostic genes were screened out (Figure 2C and Supplementary Table S3). Concurrently, 107 DEGs were identified by differential analysis between cancer and normal prostate tissue (Figure 2D and Supplementary Table S3), and a total of 44 candidate “E2F target”-related prognostic genes were obtained by intersecting these DEGs with the previously screened genes (Figure 2E and Supplementary Table S3). Subsequently, the 44 candidate genes were further screened by a LASSO Cox regression model, optimized when the minimum λ value was 0.03035, with four genes remaining, including MAX dimerization protein 3 (MXD3), polo-like kinase 1 (PLK1), EPH receptor A10 (EPHA10), and kinesin family member 4A (KIF4A) (Figures 2F,G). Figure 2H and Supplementary Table S3 show the coefficients of these four genes, and the final risk score was calculated as follows:

Risk score = MXD3* 0.579083868008915 + PLK1* 0.258020225608399 + EPHA10* 0.211234188276561 + KIF4A * 0.266044684281992.

As the risk score increased, the gene expression of the signature rose, the number of patients with BCR in the training cohort increased significantly, and the BCR-free survival time decreased (Figure 3A). The genetic correlation analysis did not show excessive high correlations (Figure 3B and Supplementary Table S4). According to the median risk score of 3.713494217, the patients were divided into high- and low-risk groups. GSEA confirmed the enrichment of the E2F target pathway genes in the high-risk group (Figure 3C). PCA and tSNE analysis showed that patients in the high- and low-risk groups could be completely distinguished from each other (Figures 3D,E). Multivariate Cox regression analysis showed that the risk score was an independent predictor of prognosis and the strongest predictor of BCR along with clinical features (HR: 2.261, 95% CI: 1.636–3.125, p < 0.001, Figure 3F). The Kaplan–Meier survival curve showed that the BCR-free survival time in the high-risk group was significantly lower than that in the low-risk group (p < 0.0001, Figure 3G). Figure 3H shows that the area under the ROC curve (AUC) of the risk score in the training group within 9 years was all higher than 0.7, which indicates that the average predictive ability of the risk score was strong. The T-ROC curve revealed that the AUC of the risk score was always higher than that of other clinical features and tended to be stable over time (Figure 3I).

The risk score of each sample was calculated in test cohorts I and II, which were divided into high- and low-risk groups according to the median risk score of the training cohort (i.e., 3.713494217). The four-gene signature showed the same expression pattern as the training cohort, and the distribution of deaths in test cohort I and metastatic patients in test cohort II increased with the risk score (Figures 4A–C). In addition, in the high-risk groups of the two cohorts, the E2F target pathway was shown to be activated (Figure 4D). In addition, the composition of the gene signature was verified to distinguish between patients at different risks and to predict poor prognosis within 10 years in the two cohorts (Figures 4E,F). The detailed risk information and gene expression of the training and test cohorts are represented in Supplementary Table S5, and the relevant information generated in the process of GSEA is shown in Supplementary Table S6.

Pearson correlation analysis was used to study the distribution of patients among different groups with various clinical feature stratifications. As shown in Figure 5A, in the training cohort, the high-risk group was closely associated with higher PSA, Gleason score (GS), and T stages. Similarly, in test cohorts I and II, the PSA, GS, and T stages tended to increase in the high-risk group (Figures 5B,C). Unfortunately, this trend was not statistically significant. In this regard, we gave a cautious explanation. In test cohort I, as OS was the primary endpoint, a shortage of other endpoints arose, which thus affected the statistically significant results. In contrast, test cohort II received radiotherapy plus androgen deprivation therapy, which, from the perspective of the E2F target pathway, may have been more beneficial to patients, thus affecting OS and the statistically significant results.

The DEGs in the high- and low-risk groups were calculated using the “limma” algorithm. Figure 6A shows all the up- and downregulated DEGs screened by using p < 0.05 as a threshold. The functional enrichment analysis of DEGs was performed, and the results of KEGG analysis showed that DNA replication, base excision repair, and mismatch repair were active in the high-risk group. Surprisingly, histidine metabolism was suppressed in this group (Figure 6B). As shown in Figures 6C,D and Supplementary Table S7, the GO analysis showed that in the high-risk group, the activation of the E2F target pathway participated in the binding of many factors and proteins and may regulate cell growth, transcription, apoptosis, and other metabolic activities by interfering with binding processes. The study also found the regulation of the Wnt and p53 signaling pathways in the high-risk group, suggesting that these pathways may be abnormally activated. In addition, the study found that the activation of the E2F target pathway was also related to remodeling of the extracellular matrix, which may enhance the invasive ability of tumor cells, promote cancer metastasis, and lead to poor prognosis. Moreover, the correlation between the risk scores and the characteristics of tumor stem cells was studied. A significant positive correlation was found between the risk score and stem cell characteristics (Figures 6E,F). All relevant information to support these results is presented in Supplementary Table S7.

Figure 7A and Supplementary Table S8 show the immune infiltration of high- and low-risk groups using different algorithms. From the ssGSEA and CIBERSORT thermograms, the immune infiltration of patients in the high-risk group was more complex with both immune activation and immunosuppression and varying degrees. However, it is worth noting that both the algorithms showed significant infiltration of T regulatory cells in the high-risk group (Figure 7B). Concomitantly, negative immune regulatory genes, such as EZH2, HAVCR1/2, and DNMT1 were also observed to be significantly expressed in the high-risk groups. When focusing on the immune subtype, it was found that the C3 type with good prognosis was concentrated in patients with low-risk scores, while the higher the risk score, the more the C1, C2, and C4 types were linked to poor prognosis (Figures 7C,D). Continuously, the IPS algorithm was used to calculate the immunophenotype of the samples with the top 10 and the last 10 risk scores, which once again verified the high expression of immune suppressor cells (SCs) and a decrease in the IPS in the high-risk group (Figure 7E and Supplementary Figure S1).

GSEA and sensitivity analysis of commonly used drugs reflected the resistance of high expression of different genes to treatment (Figures 8A,B and Supplementary Table S9). MXD3 had certain sensitivity to commonly used drugs, although not that strong. In addition, sensitive drugs were screened using the CellMiner database (Supplementary Table S9), and the results showed that -tinib drugs were more sensitive to MXD3 and EPHA10, suggesting that they may become potential new therapeutic targets.