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

Abstract

Background: 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.

Methods: 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.

Results: 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.

Conclusion: 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

Study Design

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.

Differential Expression analysis

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.

GEPIA Survival analysis

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.

Gene Ontology (GO) and Pathway analysis

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.

c-BioPortal analysis

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.

Copy Number analysis

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.

Immune Infiltration analysis

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.

Molecular Docking

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.

Molecular Dynamics Simulation

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

Identification of DEGs

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).

Table 1. Comparison of Upregulated and Downregulated Genes Across Prostate Cancer Groups.

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

Survival analysis of the differentiated genes

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).

Table 2. DEGs Associated with OS and RFS Across Prostate Cancer Groups.

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.

Figure 1. Heatmap for DEGs Associated with both OS and RFS.

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

Pathway analysis of significant genes

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).

Table 3. Go Function and Pathway Enrichment Analysis of Common Genes.

	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).

Table 4. KEGG Pathways Analysis of 7 Common Genes in Prostate Cancer.

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

Validation of survival associated DEGs by cBioPortal analysis

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).

Table 5. Survival Analysis of the Survival Associated DEGs Identified in Prostate Cancer from cBioPortal Database.

Number of patients	Gene Symbol	Stage	Risk status** (Number of patients)	Survival/deceased status** (Number of patients)	P Value	Disease free survival/recurrence status** (number of patients)	P Value
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).

Figure 2. ROC Curve of PLK1 Expression..

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.

PLK1 expression was correlated with immune cell infiltration

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).

Figure 3. Correlation between PLK1 Expression and Immune Infiltration.

Copy Number Analysis

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).

Table 6. Distribution of PLK1 Copy Number Alterations in High- and Low-expression Groups.

		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).

Figure 4. Spearman Correlation of PLK1 Expression with Copy Number in A) T2 and B) T4 Stage.

Drug sensitivity analysis

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

Figure 5. Relationship between GDSC Drug Sensitivity and PLK1 mRNA Expression.

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

Molecular dynamics simulation

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.

Figure 6. Molecular Dynamics Simulation Demonstrating A) RMSD B) RMSF, and C) Radius of gyration of Rapamycin–PLK1 Complex.

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.