Integrated Bioinformatics Reveals PLK1 as a Driver and Rapamycin as a Drug Candidate in Stage III–IV Prostate Cancer

Prostate cancer is one of the most commonly diagnosed malignancies in men worldwide. This study aimed to identify differentially expressed genes associated with prostate cancer progression by integrating tumor stage and prostate-specific antigen (PSA) levels. Microarray data from the GEO dataset (GSE116918) were analyzed using GEO2R based on gene expression levels, tumor stage, and PSA categories. The dataset comprised 51 stage I, 76 stage II, 92 stage III, and 4 stage IV tumors. Survival-associated genes were identified using GEPIA2. Receiver operating characteristic (ROC) curve analysis was performed using OriginPro software. Copy number variation analysis and immune infiltration profiling were performed using GISTIC and TIMER tools. Molecular docking and dynamic simulations were performed using PyRx and NAMD, respectively. PLK1 expression was significantly higher in stage IV than in stage II prostate cancer, indicating its potential role in tumor progression. Survival and ROC analyses further supported the prognostic significance of PLK1 in this context, demonstrating its strong predictive accuracy in advanced-stage disease. In addition to its association with disease stage and PSA levels, PLK1 expression was linked to molecular characteristics, including CpG methylation and copy number alterations, as well as immune cell infiltration involving CD8⁺ and CD4⁺ T-cells, macrophages, dendritic cells, and monocytes. Drug screening highlighted rapamycin as a potent PLK1 inhibitor, with strong binding affinity and structural stability, as confirmed by docking and molecular dynamics simulations. These findings suggest that PLK1 may serve as a stage-specific biomarker and therapeutic target, particularly in advanced-stage prostate cancer.

 Introduction

Prostate cancer (PCa) is a major health concern worldwide. It is the second most frequently diagnosed cancer in men and ranks fourth among all cancers diagnosed worldwide [1]. Based on current data, regions such as Northern Europe are seeing the highest number of cases, with an estimated incidence rate of approximately 82.8 per 100,000 men. On the other hand, Southern Africa seems to have the highest death rate from PCa, roughly around 29.7 per 100,000 men, which is quite alarming [2]. Although the exact reason for prostate cancer is not completely understood, several risk factors have been linked to it, such as age, family background, genetic influence, and ethnicity [3]. The disease tends to develop in a stepwise manner. In the early stages (T1–T2), the tumor usually remains inside the prostate gland. However, as it worsens (T3–T4), it starts spreading outside the nearby tissues [4]. In some cases, prostate cancer progresses to metastatic castration-resistant prostate cancer (mCRPC), where the cancer spreads to other parts of the body, such as the bones or lymph nodes [5]. Identifying the genetic and molecular factors underlying these stages is important for improving diagnosis and treatment. Over time, scientists have noticed several changes at both the molecular and cellular levels, which seem to play a significant role in the development and progression of prostate cancer [6]. Identifying and examining differentially expressed genes (DEGs) at various stages of prostate cancer can help us understand how the tumor grows and how it might react to treatments. However, few studies have examined DEGs in detail across all stages of PCa. For instance, the gene p53 tends to show higher expression in patients with more advanced stages of prostate cancer and higher Gleason scores, which usually indicates a more aggressive disease [7]. Similarly, changes in the expression levels of certain miRNAs have also been associated with differentiating between benign and malignant forms of prostate cancer [8]. Further DEG analysis can ultimately lead to improved overall survival (OS) in patients with prostate cancer at different stages of the disease.

Prostate-specific antigen (PSA) is a widely used diagnostic biomarker for PCa; however, its relationship with DEGs remains poorly understood. It remains unclear whether elevated PSA levels drive differential gene expression or whether DEGs contribute to increased PSA production. For example, the B4GALNT4 gene, known to be upregulated in prostate, breast, and lung cancers, plays a role in PSA glycosylation during prostate carcinogenesis [9]. Currently, bioinformatic approaches have yielded the production of public databases and improvements in in silico tools, which have ultimately aided in the analysis of DEGs [10]. The availability of high-throughput genomic datasets, such as The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets, has also enabled the identification of DEGs. The broad use of these repositories has resulted in the enhancement of precise therapeutic targets for cancer [11].

This study aimed to investigate the prognostic relevance and molecular associations of DEGs in prostate cancer, with a particular focus on their expression across different disease stages. We further explored the immune and genomic landscape linked to these DEGs and identified potential therapeutic inhibitors, ultimately evaluating their potential as stage-specific biomarkers and therapeutic targets, especially in advanced-stage prostate cancer.

Materials and Methods <italic id="italic-1">Study Design<italic id="italic-2"/></italic>

This study utilized an in-silico approach to analyze microarray data from the Gene Expression Omnibus (GEO) dataset GSE116918. The analysis, conducted using GEO2R, aimed to identify potential biomarkers and therapeutic targets for prostate cancer in patients treated with radical radiation therapy and androgen deprivation therapy. The dataset comprised 51 stage I, 76 stage II, 92 stage III, and 4 stage IV tumors, each associated with PSA levels and varying stages of prostate cancer. Subsequently, DEGs were sorted for upregulated and downregulated genes and filtered according to their expression. Furthermore, the genes that had a significant effect on survival were sorted and validated by cBioPortal data analysis. Receiver Operating Characteristic (ROC) curve analysis was performed for the selected genes, and genes with a high Area Under the Curve (AUC) were prioritized. In addition, copy number variations were assessed and immune cell infiltration associated with the selected genes was analyzed. Molecular docking and dynamics simulations were performed to assess the interactions between the selected genes and their corresponding drug candidates.

<italic id="italic-0dbd89ab2448212db66ed3879369897a">Differential</italic> <italic id="italic-f73b472780014237144c335e23920343" /> <italic id="italic-3">Expression</italic> <italic id="italic-4" /> <italic id="italic-5">analysis<italic id="italic-6"/></italic>

To identify DEGs, the samples were categorized into four groups based on tumor stage and PSA levels: 1) Group 1: T1–T2 stage with PSA < 50 ng/ml 2) Group 2: T1–T2 stage with PSA 50 ng/ml 3) Group 3: T3–T4 stage with PSA < 50 ng/ml 4) Group 4: T3–T4 stage with PSA > 50 ng/ml. Comparative analysis between these groups was performed using the GEO2R tool. DEGs were identified based on a log fold change (logFC) threshold, with genes showing logFC > 1 classified as upregulated and those with logFC < –1 classified as downregulated [12]. Statistical significance was set at p ≤ 0.05.

<italic id="italic-7">GEPIA Survival analysis<italic id="italic-8"/></italic>

Overall survival (OS) and disease-free survival (DFS) analyses of the genes identified from the GEO dataset were conducted using the Gene Expression Profiling Interactive Analysis (GEPIA) platform [13, 14]. For each gene, log-rank p-values were calculated for both overall survival (OS) and disease-free survival (DFS). Genes with a log-rank p-value of ≤ 0.05 were considered statistically significant. These significant genes were subsequently filtered and organized for further analysis.

<italic id="italic-9">Gene Ontology (GO) and Pathway analysis<italic id="italic-10"/></italic>

Functional enrichment analyses were performed using the WEB based GEne SeT AnaLysis Toolkit (WebGestalt), which includes GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses [15]. Gene ontology analysis yielded the biological function, cellular component, and molecular function of the selected genes.

<italic id="italic-11">c-BioPortal analysis<italic id="italic-12"/></italic>

To assess the clinical relevance of common DEGs associated with both OS and DFS, gene expression data for prostate cancer were retrieved from The Cancer Genome Atlas (TCGA) using the cBioPortal platform (https://www. cbioportal.org/datasets). The dataset included American Joint Committee on Cancer (AJCC) pathological staging, survival data, patient IDs, and gene expression. Based on the median ± 2 standard error of the mean (SEM), a cutoff value was established to classify patients into high- and low-risk groups according to their PLK1 expression levels. Patients in the high-risk group, marked by elevated PLK1 expression, exhibited poorer overall survival than those in the low-risk group. Subsequently, each selected gene was analyzed using the cBioPortal platform. ROC curves for the gene expression levels were generated using OriginPro version 2024. From these ROC analyses, the AUC was calculated, and PLK1 with high AUC values and a p-value < 0.05 was selected for further investigation.

<italic id="italic-9e031c7a9413e4f20004e7534c23bbab">Copy Number analysis<italic id="italic-a1e0dfe24f2bb40d60ed61f635717267"/></italic>

To comprehensively investigate the relationship between copy number variations (CNVs) and PLK1, a CNV analysis was performed. Briefly, GISTIC was used to identify regions of genomic copy number gain or loss, from which the CNV profiles of PLK1 were extracted. CNVs were classified into five categories: −2 (homozygous deletion), −1 (heterozygous loss), 0 (diploid), 1 (single copy gain), and 2 (high-level amplification). To explore the association between CNVs and gene expression levels, Spearman’s rank correlation analysis was conducted, and correlation plots were generated using OriginPro version 2024. Additionally, Fisher’s exact test was used to assess the statistical association between CNV status and prostate cancer stage, with significance determined at a p-value threshold of ≤ 0.05.

<italic id="italic-883b91216450ae984c5487215a10186f">Immune Infiltration analysis<italic id="italic-2cb20efecc44427ae7935e90f247dfdd"/></italic>

Tumor-infiltrating immune cells play a crucial role in evaluating the effectiveness of cancer therapies and their impact on patient prognosis. The correlation between PLK1 expression and tumor-infiltrating immune cells was analyzed using the TIMER 2.0 database (http://timer. comp-genomics.org/) [16]. Heatmaps were created to visually depict the correlation between PLK1 expression and the infiltration levels of various immune cell types based on Spearman’s rho values.

<italic id="italic-ae7a7121d4abcb17377782fd825a73e0">Molecular Docking<italic id="italic-9c3b28410cecc557e520af22124b7fbf"/></italic>

Inhibitor molecules targeting PLK1 were identified using the Gene Set Cancer Analysis (GSCA) database. Drugs that showed a negative correlation with the expression of the selected genes were shortlisted, and their chemical structures were obtained from the PubChem database. The three-dimensional structure of PLK1 (PDB ID: 2RKU) was retrieved from the Protein Data Bank (https://www. rcsb. org/). The protein structure was prepared by removing non-essential residues and heteroatoms using UCSF Chimera (version 1.18). Molecular docking was performed using PyRx, a virtual screening tool, to evaluate the interaction between the selected drugs and PLK1. The docked complex with the most negative binding affinity was selected for further molecular dynamics simulations.

<italic id="italic-60be7239cae2126a83406b7837fdbc82">Molecular Dynamics Simulation<italic id="italic-d17e3c75557d9392fd5b15a5e85d1ef2"/></italic>

A molecular dynamics simulation was conducted to evaluate the interaction between rapamycin and PLK1 using the NAMD3 simulation engine, with the CHARMM force field applied for the parameterization. The system preparation and simulation protocols were executed using built-in NAMD3 commands. Key structural parameters, including the root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (Rg), were computed and visualized using the Bio3D v2.3.0 package (http://thegrantlab.org/ bio). The simulation was performed for a total duration of 50 ns. The binding free energy (ΔG) of the drug–protein complex was estimated using the NAMD Energy plugin, following the methodologies previously reported by our group [17, 18].

Results <italic id="italic-96f6750811b2a1bf659c199b31cc6f1e">Identification of DEGs<italic id="italic-3d11b4577213feda761c866600206e3f"/></italic>

Microarray analysis was performed to identify DEGs across four distinct groups of prostate cancer, categorized based on disease stage and PSA levels. After preprocessing the data, a total of 16,365, 24,777, 24,893, and 25,253 DEGs were identified in Groups 1, 2, 3, and 4, respectively. To refine the selection of DEGs, stringent filtering criteria were applied, including a log fold change (logFC) threshold of ≥1 for upregulated genes and ≤ -1 for downregulated genes, with a false discovery rate (FDR) of <0.001. Comparative analysis between the groups revealed distinct gene expression patterns. In the Group 1 vs. Group 2 comparison, 256 DEGs were identified, of which 229 were upregulated and 37 were downregulated. Similarly, the comparison between Group 2 and Group 4 yielded 307 DEGs, including 260 downregulated and 47 upregulated genes. In contrast, the comparison between Groups 1 and 3 exhibited only six DEGs, with two upregulated and four downregulated genes. Finally, the comparison between Groups 3 and 4 revealed 13 DEGs, comprising five upregulated and eight downregulated genes. These findings highlight the significant molecular alterations associated with different stages of prostate cancer and PSA levels, suggesting potential biomarkers for disease progression and stratification (Table 1).

Groups Upregulated Downregulated
Group 1 &2 FAM110B, FAM171B, NEU, PACRG, CD209, NEDD9, YY1,A2M-AS1, VTA1, UBE2E3, CHRM3, ICA1L, NECAB1, RORB, KIAA1715, OR2L13, KRAS, MAOB, ABCC9, C3orf38, TMC1, ONECUT2, SLC38A11, GNAI1, MCMBP, NDUFAF6, PMCHL2, WIPF1, NUDT3, STX3, EMCN, BCOR, SLMAP, CD3G, PAIP1, KCNG3, GNAI1, AP4E1, TNPO1, A2M-AS1, ZNF429, GABRB3, LEPR, FAM118A, RBM33, IL24, PHF14, PNMA3, SETX, EGFR, UBE2E3, JAG1, NELFA, EIF3D, SIK1, CC2D2A, GRIK1, SULT1C2, MAOB, CXCL11, CGGBP1, GPR34, NDUFB2, MID1, PDSS2, HNRNPK, WDR89, UAP1, CAMK2N1, ZNF608, PPP2R1B, ERG, NEU3, FAM63B, CDH7, MFSD1, CEP97, NUFIP2, KHDRBS1, RPS6KA6, FAM98B, AGMAT, RORB, TIMM8B, QDPR, NUDT21, SSTR1, URB2, PDS5A, FRZB, SLC41A2, NPY, TMEFF2, PLIN1, UBE3A, ATP6V0E1, RPLP1, SSTR1, SMAD2, AHCTF1, SLC38A11, TRDC, ALOX15, FGD6, SPOCK1, PPM1H, DLX2, SAMD5, DLX1, BRD7, HLA-DRB1, FPR3, PHOSPHO2, LNP1, UBE2E3, ZNF813, ACER3, ZNF880, SLC38A11, NAPEPLD, AMH, C12orf66,ST8SIA6-AS1, UBE2H, IL13RA1, G3BP2, PTCHD1, LAMB1, RBM12B, KLF8, HIATL1, GALNT2, PDSS2, ENO1, GIMAP8, FBN1, LINC00894, LIN7C, HMGN2P46, ITGA6, ZNF610, DCUN1D1, AHCYL1, RBM8A, TOX, SNRPA1, SPINT1, HINT1, FAM135A, KPNA2, BAIAP2L2, MME, ORC4, C11orf70, FAM198B, SEC22B, GREM1, SP140L, MANEA, CHRM3-AS2, FAM161B, LINC00883, UTP3, PPA2, TPO, ACAP2, LRRN1, GPR176, MLLT10, NECAP2, MOB3B, LINS, ZNF577, HIPK2, TDRD1, CDO1, KCNS3, MYO6, NAAA, CGGBP1, PRPF8, CERS6, TTTY15, TRBC2, CHRM3, AHCTF1, CTNND2, KCNN2, SERBP1, F5,THSD7A, IFRD1, POLR2M, OSGIN2, SCUBE2, PI15, FAM134B, UCHL5, REL, DOCK7, SCGN, WASF1, CCDC113, GP2, LINC01006, CTPS1, UBE2D2, CPNE4, ZNF214, MED13, EIF3F, SAMD9, LINC00963, USP34, CASP2, MCTP1, RBAK, SRP14, PLSCR1, FLYWCH1, DHODH, AFAP1L2, TTTY14, ENAH, KCMF1, SLC39A8, SLC30A4, CCDC18, KLHL20 RASSF6, TEKT4P2, MAP4K3, C1QTNF9B-AS1, TRPM8, SCD, TMED5, TRPM8, NIPAL3, ALOX15B, FEN1, WWOX, PLEKHS1, TJP2, BANK1, MIPEP, NPM1, EIF3D, ERH, ABCC4, DHFRP1, TMED10, KARS, TNFRSF19, ANTXR2, ZNF826P, SPON2
Group 2&4 PDLIM5, TEKT4P2, NPM1, TEKT4P2, MXI1, NEAT1, EFCAB2, PARK7, DLG1, PILRB, ELL2, LRRC46, SCD, SNTB2, IK, HERC3, AMACR, MSR1, FEN1 /// FEN1P1, CANX, SLC30A7, UBE3B, C1QTNF9B-AS1, CHN2, CKMT1B, EZH2, CST1, TIRAP, ACTR3, BANK1, TRPM8, ZNF826P, BBC3, MAP4K3, PLK1, TCEB2, TRPM8, DMXL1, SART3, ITCH, MMP11, ZNF826P, HNRNPM, TMPRSS11D, UGT2B15, SLC4A4, LINC01088 CD209, SUPT3H, ZNF608, ICA1L, UBE3A, NEDD9, SPINT1, EMCN, SKAP2, LSAMP, BRD7, PACRG, BCL2L14, WDR89, RHOA, PTGER3, CLK3, GALNT2, NUFIP2, CDO1, MLLT10, USP6NL, NEU4, PCDH17, PLEKHA3, EIF3I, ITPR1, RBM33, FAM118A, ERG, IL20RA, TRDC, PLD5, ST8SIA6-AS1, FOXP2, MID1, EGFR, BRWD1, RORB, BRD7, WIPF1, ATP6V0E1, NECAB1, HOXD11, SMARCA1, TPO, MME, PIK3C2A, ABCA5, SPATA6, PPM1H, PIGB, SCD5, FGD6, THSD7A, ZNF880, ATL2, HMGN2P46, ZBTB20, MFSD1, LINC01006, APC, CCNDBP1 /// TMEM62,S CGN, HMGN2P46, NDUFAF6, PRDM8, PART1, WNT3, SLC38A11, LINC00161, JAG1, STON1-GTF2A1L, PNMA3, CD38, DPY19L1P1, SMAD2, CCPG1, PRDM8, PTCHD1, MAP7D3, FAM98B, GABRB3, SLC38A11, DCUN1D1, TMEM126A, MIPEPP3, FRK, LINC00883, BAIAP2L2, NUDT21, SLC25A28, CD302, ZNF711, SBSPON, SLC38A11, KLHL20, KHDRBS1, HMGN2P46, PI15, CA3, UAP1, FAM161B, MPP7, PLSCR1, TNPO1, SORBS1, FAM161A, SLC38A11, PAK1, RLIM, STRN, RNASE1, FAM149A, HOXD1, DIS3L2, STK33, ZNF765, TUBBP5, MME, PYROXD2, GIMAP8, BCOR, PLCH2, ABCC9, TMEFF2, ARL17B, SF3A1, SLC30A4, FGFR2, LINS, EPHA7, POU5F1P3, NPY, ANO5, FUT9, TMEM59, KMT2A, RNLS, CHRM3-AS2, KHDRBS1, AFAP1L2, IFIH1, ZNF214, ORC4, RBM12B, VSTM4, EPHA4, QDPR, PLG, LPAL2, SLC22A3, ZNF587B, MCMBP, FPR3, LSAMP, IGF1, APOL3, NR3C1, FAM161B, TMC1, TMEFF2, SLMAP, MBD4, PHF14, RBM8A, SP140L, CD59, KLHL20, PHKB, NFE2L2, AFAP1L2, TYW5, NEBL, RGS17, MTRF1L, UBE2FP1, AGMAT, EPS15L1, AMH, ZBTB20, ZNF404, MKL2, NFXL1, CHRM3, EPB41L3, ACER3, FAM89A, BNC2, REL, SLC11A2, HLA-DRB1, HLA-DMB, ZNF813, OSGIN2, IL13RA1, FGF13, ARL17B, KCNG3, ZNF776, SBSPON, UBE3A, SEMA6A, EXOSC6, SOX5, PIP4K2B, SOX4, HGD, PLSCR1, IFRD1, C3orf38, SH3PXD2A, CLUL1, DLX2, PTGIS, SNRPA1, KLF8, BCL2L14, PDE8B, HIVEP3, ACVR2B-AS1, TLR4, RPL24, ONECUT2, MAOB, PLG, GFRA1, WRN, PPP2R1B, ATL2, LINC00894, MLH3, CACNA2D1, CHRM3, MED13, ITPR1, RBMS2, GRIK1, MTFMT, CELF2, FAM135A, SLC38A11, DOCK7, RBAK, ZNF557, EIF4G2, HDAC1, TPRXL, PAGE4, FBN1, TMEFF2, AFAP1L2, HNRNPK, PREX2, C11orf70, ACPP, SNAPC5, C6orf203, KIF21A, EFCAB1
Group 1&3 VCAN, MATN3 MSMB, ANPEP, TPRXL /// CD81-AS1 /// NIFK-AS1, LTF
Group 3 &4 LRFN1, AMACR, DCTN4, TFAP2C, CEACAM6 ITPR1, BCAS1, KRT15, PIGR, NPY, ERG, SCGN, HLA-C

<italic id="italic-adb7c9b1268ca8794c9ed6ea0279e820">Survival analysis of the differentiated genes<italic id="italic-885c3d77886d1b9300fcb4b7b79228fe"/></italic>

Analysis of DEGs associated with OS and recurrence- free survival (RFS) across different prostate cancer groups, categorized by disease stage and PSA level, revealed distinct expression patterns. In the comparison between Groups 1 and 2, nine genes were upregulated (UOR) and one gene was downregulated (DOR) with respect to OS, whereas 30 genes were upregulated and four genes were downregulated in association with RFS. Notably, only two genes were commonly upregulated and one gene was downregulated across both OS and RFS. For Group 2 vs. Group 4, only two genes were upregulated and 12 were downregulated in association with OS, whereas 10 genes were upregulated and 45 were downregulated in relation to RFS. Among these, one gene was commonly upregulated and three genes were downregulated across both OS and RFS. A comparison between Groups 1 and 3 showed no DEGs associated with OS, and only one gene was downregulated in relation to RFS. No common DEGs were observed between the OS and RFS groups. Similarly, in the Group 3 vs. Group 4 comparison, no genes were differentially expressed in association with OS, whereas one gene was upregulated and one gene was downregulated in relation to RFS. No overlapping genes were identified between OS and RFS in this group (Table 2).

Groups Number of genes (OS)   Number of genes (RFS)   Number of genes (OS & RFS)
  UOR DOR UOR DOR UOR DOR
Group 1&2 MCMBP, NELFA, GRIK1, AMH, BAIAP2L2, UTP3, IFRD1, POLR2M, EIF3F ABCC4 A2M-AS1, CHRM3, NECAB1, GNAI1, NDUFAF6, GNAI1, A2M-AS1, FAM118A, JAG1, NELFA, NDUFB2, SSTR1, RPLP1, SSTR1, PPM1H, FPR3, ZNF813, ZNF880, LAMB1, LIN7C, ZNF610, RBM8A,SNRPA1, KPNA2, NAAA, CHRM3, IFRD1, LINC01006, FLYWCH1, CCDC18 FEN1, ABCC4, TMED10, KARS NELFA, IFRD1 ABCC4
Group 2&4 CKMT1B, PLK1 BAIAP2L2, FAM161A, ZNF587B, MCMBP, MBD4, AMH, IFRD1, CLUL1, GRIK1, ZNF557, PAGE4, KIF21A NEAT1, PILRB, MSR1, FEN1 /// FEN1P1, EZH2, BBC3, PLK1, MMP11, HNRNPM, UGT2B15, BCL2L14, PTGER3, FAM118A, NECAB1, PPM1H, ZNF880, ZBTB20, LINC01006, NDUFAF6, PART1, JAG1, CD38, CCPG1, MIPEPP3, SBSPON, MPP7, SORBS1, RNASE1, DIS3L2, POU5F1P3, ANO5, RNLS, PLG, FPR3, IGF1, MBD4, RBM8A,CD59, NFE2L2, EPS15L1, ZBTB20, CHRM3, ZNF813, SBSPON, PIP4K2B, IFRD1, SNRPA1, BCL2L14, PLG, GFRA1, CACNA2D1, CHRM3, CELF2, PAGE4, ACPP PLK1 MBD4, PAGE4, IFRD1
Group 1&3 0 0 0 MSMB 0 0
Group 3&4 0 0 LRFN1 BCAS1 0 0

A heatmap of the genes affecting both OS and RFS is shown in Figure 1.

These findings highlight significant differences in gene expression profiles across different prostate cancer stages and PSA levels, particularly concerning survival and recurrence outcomes.

<italic id="italic-ea029935a567e90524a12885ac46f7d0">Pathway analysis of significant genes<italic id="italic-41ed2028c765298877123d8678128d79"/></italic>

GO functional and KEGG pathway enrichment analyses were performed to identify potential target genes. In the biological process category, the enriched GO terms for the target genes included regulation of anaphase-promoting complex-dependent catabolic process, myoblast fate commitment, female meiosis chromosome segregation, mitotic nuclear membrane disassembly, double-strand break repair via alternative non-homologous end joining, depyrimidination, regulation of organelle assembly, positive regulation of ubiquitin-protein ligase activity, regulation of protein localization to the cell cortex, and nuclear membrane disassembly. In the cellular component category, the enriched terms included outer kinetochore, spindle midzone, mitotic spindle pole, synaptonemal complex, synaptonemal structure, spindle microtubule, condensed nuclear chromosome, sarcoplasm, postsynaptic density membrane, and postsynaptic specialization membrane. In the molecular function category, the enriched terms included anaphase-promoting complex binding, DNA N-glycosylase activity, hydrolase activity (hydrolyzing N-glycosyl compounds), DNA endonuclease activity, DNA nuclease activity, endonuclease activity, hydrolase activity (acting on glycosyl bonds), nuclease activity, magnesium ion binding, and catalytic activity on DNA (Table 3).

  Biological Process Biological Process          
  Method-1 Method-1          
Gene Set Description   Size Expect Ratio P Value FDR
GO:1905784 regulation of anaphase-promoting complex-dependent catabolic process   5 0.001786 559.9 0.001785 1
GO:0048625 myoblast fate commitment   6 0.0021432 466.58 0.0021416 1
GO:0016321 female meiosis chromosome segregation   6 0.0021432 466.58 0.0021416 1
GO:0007077 mitotic nuclear membrane disassembly   7 0.0025004 399.93 0.0024982 1
GO:0097681 double-strand break repair via alternative nonhomologous end joining   7 0.0025004 399.93 0.0024982 1
GO:0045008 depyrimidination   8 0.0028577 349.94 0.0028547 1
GO:1902115 regulation of organelle assembly   244 0.087158 22.947 0.0030333 1
GO:1904668 positive regulation of ubiquitin protein ligase activity   9 0.0032149 311.06 0.003211 1
GO:1904776 regulation of protein localization to cell cortex   9 0.0032149 311.06 0.003211 1
GO:0051081 nuclear membrane disassembly   11 0.0039293 254.5 0.0039234 1
    Cellular Component          
Gene Set Description   Size Expect Ratio P Value FDR
GO:0000940 outer kinetochore   20 0.0067571 147.99 0.0067391 1
GO:0051233 spindle midzone   37 0.012501 79.995 0.012438 1
GO:0097431 mitotic spindle pole   38 0.012839 77.89 0.012772 1
GO:0000795 synaptonemal complex   40 0.013514 73.996 0.01344 1
GO:0099086 synaptonemal structure   40 0.013514 73.996 0.01344 1
GO:0005876 spindle microtubule   79 0.026691 37.466 0.026399 1
GO:0000794 condensed nuclear chromosome   81 0.027366 36.541 0.02706 1
GO:0016528 sarcoplasm   93 0.031421 31.826 0.031017 1
GO:0098839 postsynaptic density membrane   99 0.033448 29.897 0.03299 1
GO:0099634 postsynaptic specialization membrane   120 0.040543 24.665 0.03987 1
    Molecular Function          
Gene Set Description   Size Expect Ratio P Value FDR
GO:0010997 anaphase-promoting complex binding   9 0.0015553 642.96 0.0015546 1
GO:0019104 DNA N-glycosylase activity   14 0.0024194 413.33 0.0024175 1
GO:0016799 hydrolase activity, hydrolyzing N-glycosyl compounds   40 0.0069124 144.67 0.0068969 1
GO:0004520 DNA endonuclease activity   41 0.0070853 141.14 0.0070689 1
GO:0004536 DNA nuclease activity   64 0.01106 90.417 0.01102 1
GO:0004519 endonuclease activity   122 0.021083 47.432 0.020936 1
GO:0016798 hydrolase activity, acting on glycosyl bonds   133 0.022984 43.509 0.02281 1
GO:0004518 nuclease activity   214 0.036982 27.04 0.03653 1
GO:0000287 magnesium ion binding   227 0.039228 25.492 0.03872 1
GO:0140097 catalytic activity, acting on DNA   241 0.041647 24.011 0.041074 1

The enriched KEGG pathways for the target genes included base excision repair, progesterone- mediated oocyte maturation, FOXO signaling pathway, oocyte meiosis, and cell cycle (Table 4).

Gene Set Description Size P Value FDR
hsa03410 Base excision repair 44 0.010575 1
hsa04914 Progesterone-mediated oocyte maturation 102 0.024429 1
hsa04068 FoxO signaling pathway 131 0.031319 1
hsa04114 Oocyte meiosis 131 0.031319 1
hsa04110 Cell cycle 157 0.037476 1
<sec id="heading-a2c94186a0bfdff533fc5cb71b4d4f7b"> <title> <italic id="italic-b4ae29482e332e52de8fa15188ffdd64">Validation of survival associated DEGs by cBioPortal <italic id="italic-760271e14035e1fc07ccebb1a3eb2ff0"/></italic> <italic id="italic-04c5ec29a920fd69635aea5b3bf02d85">analysis<italic id="italic-9aef8b034bf9519bae5564dfeba07358"/></italic>

The prostate cancer dataset was retrieved from cBioPortal, and RSEM data for genes influencing both OS and RFS were extracted from the dataset. Owing to the unavailability of PSA levels in the dataset, patients were stratified into two groups based on tumor stage: Group A (T1/T2) and Group B (T3/T4). Further classification into high- and low-risk groups was performed using the cut-off value (median ± SEM) of gene expression, considering their impact on OS and RFS.

None of the genes were significantly associated with the different stages of Pca (Table 5).

Number of patients GeneSymbol Stage Risk status**(Number of patients) Survival/deceased status**(Number of patients) PValue Disease free survival/recurrence status** (number of patients) PValue
195 PLK1 T2 High risk Living (71) 0.93 Censored (66) 0.39
    CO: 53.23±41.3   Deceased (1)   Progressed (6)  
      Low risk Living (80)   Censored (77)  
        Deceased (1)   Progressed (4)  
    T4 CO: High risk Living (3) 1 censored (1) 1
        Deceased (0)   Progressed (2)  
      Low risk Living (0)   Censored (0)  
        Deceased (0)   Progressed (0)  
195 MBD4 T2 High risk Living (77) 0.97 Censored (71) 0.94
        Deceased (1)   Progressed (7)  
      Low risk Living (80)   Censored (74)  
        Deceased (1)   Progressed (7)  
    T4 High risk Living (2) 1 Censored (2) 0.4
        Deceased (0)   Progressed (0)  
      Low risk Living (3)   Censored (1)  
        Deceased (0)   Progressed (2)  
195 IFRD1 T2 High risk Living (80) 0.99 Censored (71) 0.09
        Deceased (1)   Progressed (10)  
      Low risk Living (79)   Censored (76)  
        Deceased (1)   Progressed (4)  
    T4 High risk Living (3) 1 Censored (1) 0.4
        Deceased (0)   Progressed (2)  
      Low risk Living (3)   Censored (3)  
        Deceased (0)   Progressed (0)  
195 PAGE4 T2 High risk Living (71) 0.85 Censored (69) 0.07
        Deceased (1)   Progressed (3)  
      Low risk Living (61)   Censored (53)  
        Deceased (0)   Progressed (8)  
    T4 High risk Living (3) 1 Censored (3) 1
        Deceased (0)   Progressed (0)  
      Low risk Living (0)   Censored (0)  
        Deceased (0)   Progressed (0)  

*OS = Overall survival; DFS = Disease-free survival; CO = Cut-off value; *: p ≤ 0.05; **Number in brackets represent number of patients

Additionally, ROC curve analysis was performed for all genes using RSEM data. Among the analyzed genes, PLK1 exhibited the highest AUC of approximately 0.71, with a statistically significant p-value of 0.025 (Figure 2).

Although no statistically significant difference in PLK1 expression was observed between the high- and low-risk groups, ROC analysis yielded an AUC of 0.71, indicating moderate predictive accuracy and supporting the potential of PLK1 as a biomarker for prostate cancer risk stratification. Based on these findings, only PLK1 was selected for further investigation.

<italic id="italic-15079860d8d79ed5f8af2f683a241168">PLK1 expression was correlated with immune cell infiltration<italic id="italic-6a5ed1118e449f5959a3cbfbf6672c9a"/></italic>

TIMER analysis revealed a positive correlation between PLK1 expression and immune cell infiltration in prostate cancer (Figure 3). Additionally, PLK1 expression was positively associated with both immune and stromal scores. Among the various immune cell types, PLK1 expression demonstrated a strong correlation, particularly with mast cells and neutrophils, suggesting its potential involvement in modulating the tumor immune microenvironment (Figure 3).

<italic id="italic-eccae8dc68fd10ddc6713887c1b02229">Copy Number Analysis<italic id="italic-696a1095c00c3bba44b516cb1d1938d2"/></italic>

The relationship between PLK1 copy number variation (CNV) and its gene expression in prostate cancer has not been fully elucidated. In the T2 stage, both high and low PLK1 expression groups exhibited copy number deletions, no change, and amplification, with the majority showing no alteration. In contrast, in the T4 stage, no deletions were observed; only a few cases in the high-expression group showed amplification or no change, whereas the low- expression group showed no CNV events. Despite these differences, statistical analysis revealed no significant association between PLK1 expression levels and copy number variation across the T2 and T4 stages (Table 6) (p = 0.27).

    PLK1-T2 stage    
  Deletion No change Amplification Total
High 4 65 2 71
Low 4 76 1 81
    PLK1-T4 stage    
High 0 2 1 3
Low 0 0 0 0

Spearman’s correlation analysis was performed using OriginPro version 2024 to assess the relationship between PLK1 expression and copy number alterations. In the T2 stage, the correlation was weak (ρ = 0.03194) with a p-value of 0.66685, indicating no significant association between the two. In the T4 stage, a moderate positive correlation was observed (ρ = 0.40618); however, the p-value of 0.24413 suggested that this correlation was not statistically significant (Figure 4A and 4B).

<italic id="italic-b8c1831cb475eb18f3cb83c0fa60ef76">Drug sensitivity analysis<italic id="italic-70fad9bc7d98b5a58a7b2a62c91f2c1d"/></italic>

From the GSCA database, 126 drugs were selected based on their negative IC50 values relative to PLK1 expression (Figure 5).

Docking analysis revealed the highest binding affinity between PLK1-rapamycin (Supplementary Table 1).

<italic id="italic-26be0f30feb5dc6677324a75522ed36d">Molecular dynamics simulation<italic id="italic-83cf1d4da3d68462f17da06e886e0d1b"/></italic>

The amino acids that formed polar interactions with rapamycin were Lys-82, Arg-57, Arg-136, Leu-132, and Leu-59. The stability and compactness of the PLK1-rapamycin complex were assessed using RMSD and radius of gyration (Rg), as shown in Figure 6. The RMSD value did not show fluctuations throughout the 50 ns simulations, with an average of 3.73 ± 0.041 nm (Figure 6A), indicating structural stability.

In addition, the root mean square fluctuations (RMSF) were calculated over 50 ns of simulation to analyze residue-level fluctuations (Figure 6B). The complex RMSF profile indicated that Lys-82, Arg-57, Arg-136, Leu-132, and Leu-59 amino acid residues participating in polar interactions with rapamycin exhibited fluctuations below 0.9 (average:0.92 ± 0.14). The overall RMSF remained stable, averaging 1.16 ± 0.1643 nm, throughout the simulation. The average radius of gyration was 4.232 ± 0.03 nm, suggesting the stability of the 3D protein structure during molecular dynamics (MD) simulations (Figure 6C). The binding energy calculations showed negative electrostatic (-6679.42 ± 5.8 kcal/mol) and van der Waals energies (-1158.11 ± 1.32 kcal/mol), resulting in a total binding energy of -6.6 kcal/mol. These findings support the high binding affinity of the inhibitor for PLK1 during the 50 ns simulation, confirming its structural stability and potential efficacy.

 Discussion

The discovery of DEGs is essential in advancing prostate cancer research, as it may lead to improved methods for early detection, prediction of disease course, and tailoring of treatment. Although numerous genes have been found to play a role in the development of prostate cancer, their value in clinical settings becomes most evident when they are associated with patient survival [8, 19]. Focusing on genes that influence both OS and RFS offers valuable opportunities for developing targeted treatments and personalized approaches in prostate cancer care. As gene expression can differ across various tumor stages, as shown in other cancer types, it is important to examine DEGs in a stage-specific manner. However, there is currently a lack of thorough analysis of DEGs across different stages of prostate cancer, representing a significant gap in the existing research. Filling this gap by identifying stage-specific DEGs is vital for enhancing the accuracy of prognosis, improving treatment precision, discovering reliable biomarkers, and gaining deeper insight into how the disease progresses, ultimately supporting more effective and individualized patient management.

In this study, we addressed this gap by examining DEGs across prostate cancer stages, revealing distinct molecular profiles linked to disease progression. By stratifying patients based on stage and PSA levels, we uncovered specific gene expression patterns relevant to prognosis. Among the survival-associated genes identified, NELFA, IFRD1, and ABCC4 were significant in early stage disease, whereas PLK1, MBD4, IFRD1, and PAGE4 were significant in advanced stages. Among the survival-associated genes identified, PLK1 emerged as a particularly notable candidate, showing strong prognostic relevance, specifically in the comparison between stages II and IV. In other stage comparisons, however, PLK1 did not show a similarly strong performance, suggesting that its role as a biomarker may be limited to specific stages of disease progression. This underscores the importance of conducting stage-wise evaluations, as markers such as PLK1 may not exhibit universal applicability throughout all phases of prostate cancer. The link between PLK1 expression and the presence of mast cells and neutrophils also suggests that it could influence the tumor immune environment, especially in more advanced stages. Furthermore, interaction studies with rapamycin showed a stable binding configuration, supported by molecular dynamics simulations and a negative binding energy of -6.6 kcal/mol. These results indicate that PLK1 could serve as both a selective prognostic indicator and a candidate for targeted therapy in prostate cancer, particularly during the transition to more aggressive stages.

Previous studies have shown that PLK1 acts as a potential prognostic marker in prostate cancer, mainly through its influence on the androgen receptor (AR) signaling pathway [20]. Blocking PLK1 interferes with mitosis and triggers cell death in multiple prostate cancer cell lines, and interestingly, this effect occurs even without relying on AR signaling [21]. Some studies have even pointed out that using PLK1 inhibitors along with abiraterone can work together to disturb mitotic processes and lead to necroptosis in cases of CRPC [22]. This suggests that PLK1 could be a valuable treatment target, especially for more advanced or AR-independent tumors. Our study supports this view but further refines it by demonstrating that the prognostic relevance of PLK1 is more evident in later stages of the disease, especially in stage II versus stage IV. Importantly, PLK1 did not show high predictive accuracy in other stage comparisons, suggesting that its clinical utility may be limited to specific phases of tumor progression. We also observed no significant correlation between PLK1 expression and copy number alterations, indicating that its role in disease progression may involve transcriptional or post-transcriptional mechanisms rather than genomic amplification. Unlike many studies that rely primarily on computational models, such as Lasso regression, to identify key genes associated with prostate cancer or CRPC, our approach integrated multiple layers of analysis [23]. By incorporating survival data, functional annotations, immune profiling, and therapeutic docking, we conducted a multifaceted assessment of PLK1’s role in prostate cancer biology. This allowed for a more comprehensive evaluation of its potential as a biomarker and drug target.

Furthermore, our findings align with existing research demonstrating the inhibitory effects of rapamycin on prostate cancer. Previous studies have shown that rapamycin suppresses prostate cancer cell proliferation by modulating cyclin D1 and potentiates the cytotoxic effects of cisplatin [24]. Given that the PI3K-AKT-mTOR pathway plays a central role in prostate cancer progression and resistance to treatments such as enzalutamide, targeting this pathway is of particular interest [25]. PLK1 has been implicated in the promotion of this signaling cascade, leading to increased androgen biosynthesis and contributing to CRPC development [25, 26]. In our analysis, molecular docking revealed a stable interaction between PLK1 and rapamycin, which was supported by simulation data and favorable binding energy. These findings suggest that targeting PLK1, particularly through modulation of the PI3K-AKT-mTOR axis, could be a promising therapeutic strategy for advanced prostate cancer.

Limitations of the study

Despite these promising findings, some limitations of this study must be acknowledged. The sample size in the GEO dataset was relatively small, which may affect the generalizability of our results. Additionally, while in silico analyses provide valuable insights, experimental validation through in vitro and in vivo studies is essential to confirm the functional relevance of the identified genes and their therapeutic potentials. Future studies should focus on validating these findings in larger, independent cohorts and exploring the mechanistic role of PLK1 in prostate cancer.

Acknowledgements

Not applicable.

Funding

Intramural funding from Parul University (RDC/ IMSL/213) was received for this study.

Author information

DB: Conceptualization, methodology, writing original draft preparation; CJ: Data analysis; MP and SP: making tables, figures, and literature search; RG: Supervision, editing the manuscript. All authors have read and approved the manuscript.

Availability of data and materials

All Datasets included in this study can be downloaded from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this study.

References 3 05 229-263 74 2024 10.3322/caac.21834 Bray F Laversanne M Sung H Ferlay J Siegel RL Soerjomataram I Jemal A CA: A Cancer Journal for Clinicians Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries 3 03 302-313 87 2025 10.1016/j.eururo.2024.11.013 Schafer EJ Laversanne M Sung H Soerjomataram I Briganti A Dahut W Bray F Jemal A European Urology Recent Patterns and Trends in Global Prostate Cancer Incidence and Mortality: An Update 2021 Ng KL Prostate Cancer The Etiology of Prostate Cancer. In: Bott SRJ, Ng KL, eds 4 10 830-844 189 2007 10.2214/AJR.07.2011 Kundra V Silverman PM Matin SF Choi H American Journal of Roentgenology Imaging in Oncology from The University of Texas M. D. Anderson Cancer Center: Diagnosis, Staging, and Surveillance of Prostate Cancer 07 3 03 1400-1415 31 2024 10.3390/curroncol31030106 Chi KN Yip SM Bauman G Probst S Emmenegger U Kollmannsberger CK Martineau P et al Current Oncology 177Lu-PSMA-617 in Metastatic Castration-Resistant Prostate Cancer: A Review of the Evidence and Implications for Canadian Clinical Practice 30 3 07 82 6 2019 10.3390/medicines6030082 Testa U Castelli G Pelosi E Medicines Cellular and Molecular Mechanisms Underlying Prostate Cancer Development: Therapeutic Implications 18 12 12 e0145322 10 2015 10.1371/journal.pone.0145322 Myers JS Von Lersner AK Robbins CJ Sang QA PLOS ONE Differentially Expressed Genes and Signature Pathways of Human Prostate Cancer 08 4 04 655 13 2022 10.3390/genes13040655 Khan MM Mohsen MT Malik MZ Bagabir SA Alkhanani MF Haque S Serajuddin M Bharadwaj M Genes Identification of Potential Key Genes in Prostate Cancer with Gene Expression, Pivotal Pathways and Regulatory Networks Analysis Using Integrated Bioinformatics Methods 09 08 730 10 2019 10.3389/fgene.2019.00730 Pudova EA Lukyanova EN Nyushko KM Mikhaylenko DS Zaretsky AR Snezhkina AV Savvateeva MV et al Frontiers in Genetics Differentially Expressed Genes Associated With Prognosis in Locally Advanced Lymph Node-Negative Prostate Cancer 02 1 07 9912 8 2018 10.1038/s41598-018-28299-9 Han J Chen M Wang Y Gong B Zhuang T Liang L Qiao H Scientific Reports Identification of Biomarkers Based on Differentially Expressed Genes in Papillary Thyroid Carcinoma 881246 12 2022 10.3389/fonc.2022.881246 Mohammad T Singh P Jairajpuri DS Al-Keridis LA Alshammari N Adnan M Dohare R Hassan MI Frontiers in Oncology Differential Gene Expression and Weighted Correlation Network Dynamics in High-Throughput Datasets of Prostate Cancer 320 3 327 8 2021 Patel MS Pratik; Gandupalli, Lithip; Gupta, Reeshu. iomedical and Biotechnology Research Journal Identification and Evaluation of Survival-associated Common Chemoresistant Genes in Cancer 02 W1 07 W556-W560 47 2019 10.1093/nar/gkz430 Tang Z Kang B Li C Chen T Zhang Z Nucleic Acids Research GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis 03 W1 07 W98-W102 45 2017 10.1093/nar/gkx247 Tang Z Li C Kang B Gao G Li C Zhang Z Nucleic Acids Research GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses 02 W1 07 W199-W205 47 2019 10.1093/nar/gkz401 Liao Y Wang J Jaehnig EJ Shi Z Zhang B Nucleic Acids Research WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs 10 10 122861 270 2020 10.1016/j.jclepro.2020.122861 Dhir A Talwar S Kaur P Malibari A Journal of Cleaner Production Food waste in hospitality and food services: A systematic literature review and framework development approach 268 53 6 2024 Spandidos Publications style Singh P PM Bhowmik D Kumari N Prajapati SK Gupta R World Acad Sci J Identification of common biomarkers affecting patient survival in cancers 63 6 2024 Prajapati S. PB Patel M Gupta R World Acad Sci J Acidic pH at physiological salinity enhances the antitumor efficacy of lenvatinib, a drug targeting vascular endothelial growth factor receptors 1168772 13 2023 10.3389/fonc.2023.1168772 Zhu H. Lin Q. Gao X. Huang X. Front Oncol Identification of the hub genes associated with prostate cancer tumorigenesis 6 41 2023 CW SN Liu AX Zhang Z Colman H Lam ET Carpten JD Singer EA et al Journal of Clinical Oncology Polo-like kinase-1 as a potential prognostic marker of prostate cancer utilizing ORIEN data 24 2973–83 32 2013 10.1038/onc.2012.309 Deeraksa A. Pan J. Sha Y. Oncogene. Jun Plk1 is upregulated in androgen-insensitive prostate cancer cells and its inhibition leads to necroptosis 2 219–238 83 2023 10.1158/0008-5472.CAN-22-1533 Patterson J.C. Varkaris A. Croucher P.J.P. Cancer Res. Jan Plk1 Inhibitors and Abiraterone Synergistically Disrupt Mitosis and Kill Cancer Cells of Disparate Origin Independently of Androgen Receptor Signaling 1 15 2025 10.1038/s41598-025-90444-y Fan C. Huang Z. Xu H. Sci Rep. Feb Machine learning-based identification of co-expressed genes in prostate cancer and CRPC and construction of prognostic models 8 1772-1784 6 2016 Imrali A Mao X Yeste-Velasco M Shamash J Lu Y American Journal of Cancer Research Rapamycin inhibits prostate cancer cell growth through cyclin D1 and enhances the cytotoxic efficacy of cisplatin 09 152258 777 2025 10.1016/j.bbrc.2025.152258 Zhang Y Li C Guo C Liu P Li C Zhao H Gong Z Biochemical and Biophysical Research Communications lncRNA HOTAIR regulates radio-resistance in squamous cell carcinoma of the tongue by Notch signaling 4 08 852-863 11 2018 10.1016/j.tranon.2018.04.011 Mao Y Xi L Li Q Wang S Cai Z Zhang X Yu C Translational Oncology Combination of PI3K/Akt Pathway Inhibition and Plk1 Depletion Can Enhance Chemosensitivity to Gemcitabine in Pancreatic Carcinoma