Definition 1. Let graph G (V, E) denote a graph with a collection of protein nodes, where $V=\{v_1,\mathrm{\dots },v_{12}\}$ and $E=\{e_{ij}\}$. G and E corresponds to the graph and edges set for protein–protein interaction (PPI) networks.

Definition 2. Given graph G (V, E), graph embedding is a mapping $f:v_i\to y_i\in {\text{𝐑}}^d\forall i\in \left[n\right]$ such that $d\ll |V|$ and the function $f$ preserve some proximity measure defined on graph G.

Each node in the graph embedding could be mapped to a low-dimensional feature vector.

Three cohorts were utilized in our study; the descriptions of these cohorts are as follows. GSE6919 and GSE30174 were obtained from the GEO database as training datasets. GSE6919 was based on the Agilent GPL92, GPL93, and GPL8300 platforms (Affymetrix Human Genome U95 Version 2 Array), submitted by Federico Alberto Monzon (2018). The GSE6919 dataset contained 504 samples, including 233 normal prostate tissues and 271 metastatic prostate tumors. GSE30174 was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), submitted by Jennifer Barb (2019). The GSE30174 dataset contained 80 samples, including 10 healthy peripheral blood and 70 non-metastatic prostate tumors. GSE16560 as a validation dataset was based on the GPL5474 platform (Human 6k Transcriptionally Informative Gene Panel for DASL), submitted by Andrea Sboner (2013), contained 281 samples including primary prostate tumors ordered by different GS. Dataset information is summarized in Table 1.

All of the training datasets were analyzed using the online tool GEO2R. Then the DEGs were calculated using the limma R package (version 3.36.5, WEHI Bioinformatics, Bundoora, Victoria, Australia) between normal and tumor samples [20]. Multiple testing was corrected by the Benjamini and Hochberg (BH) [21] method to obtain the adjusted p value. The cutoff values for screening DEGs were set at an adjusted p 0.05 and $|$log2FC$|$ $>$1.5.

GO analysis is a widely used technique for annotating genes and gene products to identify the biological process (BP), cellular component (CC), and molecular function (MF) [22]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database serves as a valuable resource for systematically analyzing gene functions and connecting genomic information with higher level functional information [23]. To ensure the success of any high-throughput gene functional analysis, it is crucial to map the user’s genes to the appropriate biological annotations within the Database for Annotation, Visualization, and Integrated Discovery (DAVID) database [24]. This serves as a foundational step for the analysis process.

To construct a graph-embedding network based on the PPI network extracted from gene expression profile data, the Search Tool for the Retrieval of Interacting Genes (STRING) database was utilized to transform gene expression data into vector input. To assess the interactive associations among DEGs, we carried out mapping of the DEGs onto STRING database, exclusively considering experimentally confirmed interactions with a combined score $\ge$0.4. Graphs are good at denoting information and relationships in PPI networks [25]. Due to the presence of identical protein complexes, proteins that possess similar topology exhibit comparable biological functions [26]. Consequently, the exploration of distinctive structural attributes among PPI networks associated with cancer might offer valuable insights into the advancement of tumors. With the aim of extracting, characterizing, and distinguishing network structures within intricate PPI networks, this research paper introduces an extensive and organized analysis of diverse graph embedding algorithms. Within the realm of graph embedding, numerous studies have been conducted, concentrating on the development of novel algorithms suitable for graphs encompassing millions of nodes and edges.

In this paper, we employed an unsupervised learning framework known as GAEs (graph autoencoders) to encode graphs or nodes into a latent vector space, and subsequently reconstructed the graph data based on the encoded information. GAEs can encode the node features and graph construction into the latent representations and decode the graph construction. So GAEs are suitable to acquire knowledge on network embedding and graph generative distributions from PPI nodes. Regarding network embedding, GAEs primarily focus on learning representations of protein nodes through the PPI graph’s structural characteristics such as the graph adjacency matrix. Regarding graph generation, diverse approaches exist, including step-by-step generation of nodes and edges or the simultaneous generation of an entire graph. Finally, the GAE results were visualized by Cytoscape software.

To assess the validity of the hub genes and drug targets identified in the training dataset, cross-validation was conducted utilizing the merged databases. The databases utilized in this study are listed in Table 2 (Ref. [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]).

Evaluation of the performance of GAEs for 10 hub genes was performed by receiver operating characteristic (ROC) analysis (Supplementary Material 1).

To further investigate the relationship between patient overall survival (OS) and the expression levels of individual hub genes, univariate Cox analysis [40, 41] was carried out to construct a prognostic model of PCa. The survival and survminer R package [42] were utilized for this analysis. The screened genes with p 0.01 in univariate Cox regression analysis were considered significant. Next, an independent investigation using the multivariate Cox proportional hazards regression analysis was performed to evaluate the influence of numerous genes as individual prognostic indicators impacting the survival of patients [43]. Finally, the stepwise method was used to select the optimized model. The risk score was calculated as follows:

$$\begin{array}{c}\hfill \mathrm{risk}\mathrm{score}=\sum \mathrm{coefficient}\mathrm{of}{\mathrm{gene}}_{\mathrm{i}}\times \mathrm{expression}\mathrm{level}\mathrm{of}{\mathrm{gene}}_{\mathrm{i}}.\hfill \end{array}$$

ROC analysis was performed to assess the prognostic risk model, and the area under the curve (AUC) was calculated. The performance of the prognostic risk model was then validated using the GEO dataset GSE16560.

After data preprocessing, 6269 of 26,696 DEGs were identified; the top 200 DEGs included 153 upregulated genes and 47 downregulated genes (Supplementary Fig. 1). The DEGs are shown in Supplementary Fig. 2.

To identify overrepresented GO categories and KEGG pathways, the obtained DEGs were subjected to analysis using the online software DAVID. The GO analysis results (Supplementary Material 2) revealed the significant enrichment of DEGs in various BPs such as signal transduction, positive regulation of transcription from RNA polymerase II promoter, and cell adhesion (Fig. 2a). Furthermore, GO analysis demonstrated the significant enrichment of DEGs in cellular compartments such as the cytoplasmic vesicle membrane, integral component of membrane, and plasma membrane (Fig. 2b). Regarding the MF, the DEGs showed enrichment in protein binding, protein homodimerization, and calcium ion binding (Fig. 2c). The overall distribution of the GO results are shown in Supplementary Fig. 3.

Fig. 2.

Significant enriched genes. (a) Biological process. (b) Cellular component. (c) Molecular function.

Table 3 presents the enriched pathways of the upregulated and downregulated DEGs. These pathways were analyzed using KEGG analysis. The upregulated DEGs demonstrated enrichment in pathways such as the cell cycle, DNA replication, progesterone-mediated oocyte maturation, p53 signaling pathway, and extracellular matrix (ECM)–receptor interaction. Conversely, the downregulated DEGs exhibited enrichment in pathways including drug metabolism, metabolism of xenobiotics by cytochrome P450, retinol metabolism, hematopoietic cell lineage, and calcium signaling. The results of KEGG enrichment analysis conducted by GSEA [44] is shown in Supplementary Fig. 4. All of the upregulated genes were significantly enriched in dilated cardiomyopathy, hypertrophic cardiomyopathy, ECM–receptor interaction, arrhythmogenic right ventricular cardiomyopathy, focal adhesion, and the transforming growth factor beta signaling pathway.

All DEGs were uploaded to the STRING database, with the analysis results shown in Supplementary Fig. 5. Based on the information available in the STRING database, a total of 6475 nodes were identified in topological form. In this paper, we present the illustrated GAEs, which utilize an encoder for extracting network embeddings and employ a decoder to ensure the preservation of nodes’ topological information through the adjacency matrix [45] as:

$${\text{PPMI}}_{v_1,v_2}=\text{max}\left(\text{log}\mathrm{}(\frac{\text{count}(v_1,v_2)\cdot \left|\text{D}\right|}{count\left(v_1\right)count\left(v_2\right)})\right),0)$$

Where $v_1,v_2\in \text{V},\left|\text{D}\right|={\sum }_{v_1,v_2}\text{count}(v_1,v_2)$ and the $\text{count}(\cdot )$ function returns the frequency that node $v$ and/or node $u$ co-occur/occur in sampled random walks.

The top 100 genes generated by GAEs are visualized in Supplementary Fig. 6. The significantly upregulated genes were ubiquitin-conjugating enzyme E2C (UBE2C), cyclin B1 (CCNB1), topoisomerase II$\alpha$ (TOP2A), Xenopus kinesin-like protein 2 targeting protein (TPX2), centromere protein M (CENPM), coagulation factor V (F5), apolipoprotein E (APOE), neuropeptide Y (NPY), and tripartite motif-containing 36 (TRIM36). The significantly downregulated genes were myosin heavy chain 11 (MYH11), filamin A (FLNA), actin alpha 2, smooth muscle (ACTA2), myosin light chain 9 (MYL9), transgelin (TAGLN), and actin gamma 2, smooth muscle (ACTG2) (Table 4).

To conduct deeper exploration of the potential hub genes, validation of these findings was carried out using the Gene Expression Profiling Interactive Analysis database. The results obtained from this analysis are illustrated in Fig. 3. Notably, the hub genes demonstrated distinct significance in both the tumor and normal groups.

Fig. 3.

Validation of hub genes by Gene Expression Profiling Interactive Analysis (red: tumor sample; grey: normal sample). PRAD, Prostate adenocarcinoma.

To confirm the hub gene expression among various cancers, the Oncomine database was employed to evaluate the expression of hub genes in tumor and normal tissues. By employing the criteria of p 0.01 and $|$log2FC$|$ $>$1.5, a total of 396, 470, 461, 426, 459, 470, 452, 399, 419, and 443 unique analyses for UBE2C, CCNB1, TOP2A, TPX2, CENPM, F5, APOE, NPY, and TRIM36, respectively, were performed in Fig. 4.

Fig. 4.

Validation of hub genes in the Oncomine database.

Some of the downregulated hub genes were enriched in ACTG2, MYH11, MYL9, FLNA, TAGLN, and ACTA2, which were validated by the Genotype-Tissue Expression database (Fig. 5). Genotype-Tissue Expression (GTEx) Program established a data resource and tissue bank to study the relationship between genetic variants and gene expression in multiple human tissues.

Fig. 5.

Significant downregulated genes in the Genotype-Tissue Expression (GTEx) database.

The Human Protein Atlas database was utilized for the validation of hub genes (Fig. 6). However, there was no pathology results for TRIM36 in this database.

Fig. 6.

Validation of the hub genes was performed using the Human Protein Atlas database.

Given the crucial role of hub genes in PCa, we analyzed the role of hub genes in predicting OS in PCa using univariate Cox regression analysis. Our findings revealed that nine genes exhibited a significant association with OS in patients with PCa (Table 5).

Finally, a stepwise multivariate Cox proportional hazards model was constructed (Table 6), and four genes were selected to build the following risk model:

$$\begin{array}{c}\hfill \text{risk score = (–0.5663 * expression level of TOP2A)}\hfill \\ \hfill \text{+ (–1.2489 * expression level of UBE2C)}\hfill \\ \hfill \text{+ (–1.4976 * expression level of MYL9)}\hfill \\ \hfill \text{+ (–0.9500 * expression level of FLNA)}\hfill \end{array}$$

Risk scores were computed for every patient in the training group, and subsequently, the patients were categorized into high- and low-risk groups. Fig. 7 illustrates the Kaplan–Meier OS curves, which clearly demonstrated that patients with high-risk scores exhibited markedly distinct survival outcomes compared to those with low-risk scores. The AUC values of the four-gene biomarker prognostic model at 1, 3, and 5 years were 0.973, 0.793, and 0.66, respectively (Fig. 8).

Fig. 7.

Kaplan–Meier plots for biomarkers in the high- and low risk groups.

Fig. 8.

Time-dependent (receiver operating characteristic, ROC) biomarker curves of 1-, 3- and 5-year overall survival.

To assess the autonomous predictive significance of a prognostic model composed of four genes, univariate and multivariate Cox regression analyses were implemented using both TCGA prostate adenocarcinoma (PRAD) cohort and GSE 16560 cohort. Univariate Cox regression analysis demonstrated that the prognostic model incorporating clinical information such as GS and pathologic stage exhibited some prognostic value (Fig. 9). The statistical significance of age and GS was nearly achieved, which prompted us to include age, GS, and the prognostic model in the multivariate Cox regression analysis (Fig. 10). The results of the multivariate Cox regression analysis showed that the prognostic model was independently associated with OS.

Fig. 9.

Forest plot of risk scores and clinical factors based on univariate Cox regression analysis.

Fig. 10.

Forest plot of risk scores and clinical factors based on multivariate Cox regression analysis.

In addition, the prognostic model’s predictive value was evaluated using the GSE16560 dataset. The dataset consisted of 280 patients, who were divided into a high-risk group (n = 190) and low-risk group (n = 90) based on the optimal cut-off value (Fig. 11). Fig. 12 illustrates the time-dependent ROC analysis results for the prognostic model’s survival prediction, which showed AUC values of 0.69, 0.58, and 0.61 at 1, 3, and 5 years, respectively.

Fig. 11.

Kaplan–Meier plots for biomarkers in the GSE16560 dataset.

Fig. 12.

Time-dependent (ROC) curve of biomarkers in the GSE16560 dataset.

All of the hub genes were uploaded to the cBioPortal to make pan-cancer analysis by TCGA PRAD dataset. The results showed that most of the hub genes had significant mutations (Supplementary Fig. 7). Furthermore, network analysis showed that the hub genes of NYP, TOP2A, and TPX2 could also serve as drug targets (Fig. 13).

Fig. 13.

Integration analysis for drug target. (a) Expression of hub genes in TCGA pan-cancer dataset. (b) Alteration frequency of hub genes in PCa. (c) Clustering of multivariate data. (d) Network of hub genes and drug targets.

In the DrugBank database, the drugs amsacrine, bicalutamide, dexrazoxane, doxorubicin, daunorubicin, enzalutamide, epirubicin, fleroxacin, mitoxantrone, teniposide, and valrubicin were PCa target drugs (Supplementary Material 3). The drug targets calculated by systemsDock showed that NY, TOP2A, and TPX2 had significant docking scores with drug targets (Fig. 14).

Fig. 14.

Molecule docking prediction for drug target genes (NY, TOP2A, and TPX2).

PCa is usually classified by the GS system. The analysis of PCa tissue is done through the GS grading system, which characterizes the tissue based on its microscopic growth pattern. This evaluation method, known as the GS system, encompasses various levels of stratification. These include GS $\le$6, 3+4, 4+3, 8, 4+5, 5+4, and 10, corresponding with Gleason Grading Group 1, 2, 3, 4, and 5, respectively [46]. Consequently, GS ranges from 4 (2+2) to 10 (5+5). A favorable prognosis can be expected for individuals with a low GS score ($\le$6), as it indicates no risk of lymphatic metastasis. Conversely, individuals with a high GS score ($>$8) are more likely to experience distal metastasis. Table 7 shows the significant deregulation of UBE2C, CCNB1, TOP2A, TPX2, CENPM, F5, APOE, NPY, and TRIM36 in PCa based on the GS. Significant gene expression of UBE2C, CCNB1, TOP2A, TPX2, CENPM, and APOE was observed in the tumor samples (Fig. 15). There was no expression of NYP, which had a GS of 10. A nomogram with clinical information indicated that the GS score showed moderate value in PCa (Supplementary Fig. 8).

Fig. 15.

Expression of hub genes according to the Gleason Score system.

PCa exhibited PSA dependence on AR signaling [47]. The relationship between PSA and the AR pathway is shown in Supplementary Fig. 9. PSA is a significant biomarker that is currently used for PCa screening and diagnosis [48]. AR splice variants may contribute to the progression of PCa. Several therapeutic drugs, including bicalutamide and enzalutamide, specifically target AR signaling in the treatment of PCa [49].