All PCa datasets were obtained from the GEO public database (http://www.ncbi.nlm.nih.gov/geo) including the GSE46602, GSE79021, GSE200879, GSE60329, and GSE71016 datasets. These datasets were selected based on their inclusion of both tumor and normal prostate tissue samples, as well as their relevance to prostate cancer biomarker discovery. These datasets were divided into a training group (GSE46602 and GSE79021) with 189 PCa cases and 63 normal control cases, and the test group (GSE200879, GSE60329, and GSE71016) with 216 PCa cases and 70 normal control cases. Detailed information about each dataset, including sample collection methods, sources, purposes, and database curators, is provided in Table 1 . Hyperlinks to the original datasets are included for further reference.

The raw data were processed into standardized data using the R package “limma”. We screened data from 189 PCa cases and 63 normal control cases in the training group and employed the “SVA” package to correct for batch effects. The threshold for DEG selection was set as |logFold Change (logFC)| > 1 and adjusted p-value < 0.05. Volcano plots for all DEGs and clustering heatmaps were generated using the R packages “ggplot2” and “pheatmap”, respectively.

The R package “WGCNA” is used to construct a gene co-expression network and identify modules related to PCa. A subset of genes with a standard deviation greater than 0.5 is selected for further analysis. The “goodSamplesGenes” function is employed to check for missing values in the data, to remove genes or samples that do not meet the quality standards, and to build a scale-free co-expression network. The “pickSoftThreshold” function is utilized to determine the optimal soft threshold. On this basis, the gene expression data matrix is converted into the corresponding adjacency matrix, which is then transformed into a topological overlap matrix to determine the gene weights and similarities. Gene modules are identified through topological overlap clustering methods. The module eigengenes (ME) are calculated and merged among similar modules, followed by the generation of a hierarchical clustering dendrogram. Subsequently, the Gene significance (GS) and Module Membership (MM) are computed within the modules through intramodular analysis.

Using the R package “clusterProfiler,” core genes were subjected to Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, with a p-value of <0.05 set as the criterion. To further explore the signaling pathways involved in core differential genes, we employed the Gene Set Variation Analysis (GSVA) method. The R packages “GSVA” and “enrichplot” were used for GSVA analysis of the differential genes, referring to the gene set file “c2.cp.kegg.Hs.symbols.gmt”.

Machine learning is a branch of artificial intelligence, and the core of machine learning is the algorithms, which can analyze and extract patterns from data and then use these patterns to make predictions or decisions without explicit programming (15). We divided the PCa dataset into a training group and a test group and included 12 machine learning methods, including SVM, Ridge, Enet, glmBoost, Lasso, plsRglm, Stepglm, RF, GBM, LDA, XGBoost, and NaiveBayes. We constructed models with 113 combinations of the 12 algorithms in the PCa training group dataset, and all constructed models were evaluated in the PCa test group dataset. For each model, we calculated its AUC values in the training and validation groups and extracted the feature genes screened by each method. Then, we ranked the predictive performance of the models based on the average AUC values to select the optimal model and the optimal model genes.

The CIBERSORT algorithm was utilized to assess the composition and abundance of 22 types of immune-infiltrating cells in both PCa and control samples. The R package “Cibersort” was employed for the immune cell infiltration analysis. Bar charts were used to visualize the proportion of each type of immune cell in different samples. Comparisons in the proportions of different types of immune cells between the PCa and control groups were visualized using the R package “ggpubr”. A heatmap depicting the correlations of the 22 infiltrating immune cells was plotted using the R package “corrplot”.

We conducted Spearman’s rank correlation analysis using R software to examine the correlation between the expression levels of the identified biomarkers and the levels of infiltrating immune cells. Finally, the results of the analysis were visualized using the R packages “ggpubr” and the “linkET” function.

We obtained expression quantitative trait locus (eQTL) data from the IEU Open GWAS project (https://gwas.mrcieu.ac.uk/), which covers 31,684 blood and peripheral blood mononuclear cell samples involving 19,942 genes. Concurrently, we obtained protein quantitative trait locus (pQTL) data from https://www.decode.com/summarydate/, encompassing 4,907 proteins in the plasma of 35,559 Icelandic individuals, and identified 18,084 sequence variations associated with plasma protein levels. We used these two datasets as exposure data for subsequent Mendelian Randomization (MR) analysis. Furthermore, the PCa outcome data “Finngen_R10_C3_PROSTATE” for MR analysis was derived from https://www.finngen.fi/en/, including 15,199 PCa cases and 131,266 control cases. In our study, single nucleotide polymorphisms (SNPs) were defined as instrumental variables (IVs), with the selection criteria for IVs being: (a) demonstrating a genome-wide significant association (P < 5 × 10^–8) (16); (b) showing an independent association (linkage disquilibrium (LD) > 10,000 kb, clump r2 < 0.001) (17). Additionally, to avoid bias caused by weak instruments, we considered those IVs with an F statistic >10 as strong instruments and reserved them for the following analysis (18). In the causal inference of Mendelian Randomization (MR), valid instrumental variables (IVs) must satisfy three fundamental assumptions: (1) the genetic variant is directly related to the exposure; (2) the genetic variant is unrelated to potential confounders between the exposure and the outcome; (3) the genetic variant does not affect the outcome through pathways other than the exposure.

This study was approved by the Ethics Committee of Wuxi No. 2 Hospital affiliated with Nanjing Medical University (2022-Y-80). The PCa specimens and their paired normal tissues were obtained from patients at Wuxi No. 2. Hospital and all participants have signed informed consent forms. All patients received no endocrine therapy before surgery, and all patients underwent radical prostatectomy. The tissues were immediately preserved in liquid nitrogen.

FastPure Cell/Tissue Total RNA Isolation Kit V2, HiScript III RT SuperMix, and ChamQ Universal SYBR qPCR Master Mix Q711-02 were all acquired from Vazyme (Nanjing, China). RIPA lysis buffer was purchased from Beyotime (Shanghai, China). Protease Inhibitor Cocktails (Cat. HY-K0010) were sourced from MedChemExpress (Shanghai, China). Anti- ARHGEF38 (Cat. ab122345) was purchased from Abcam(Shanghai, China). Antibodies against UT-B (Cat. 25962-1-AP), NEFH (Cat. 18934-1-AP), MSMB (Cat. 15888-1-AP), KRT23 (Cat. 24049-1-AP), KRT15 (Cat. 10137-1-AP) and Beta Actin (Cat. 66009-1-Ig) were purchased from Proteintech (Wuhan, China). Secondary antibodies (Cat. SA00001-1and Cat. SA00001-2) were purchased from Proteintech (Wuhan, China).

Total RNA was isolated from cells with FastPure Tissue Total RNA Isolation Kit V2 and reverse-transcribed into complementary DNAs (cDNAs) with HiScript III RT SuperMix. The cDNAs were amplified based on the standard qPCR protocol with ChamQ Universal SYBR qPCR Master Mix in an Applied Biosystems 7500 Fast Real-Time PCR System (Thermofisher, US). Quantitative PCR was conducted at 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. GAPDH was used as the internal control and the results were calculated using 2^-ΔΔCT method. The sequences of primers can be seen in Supplementary Table 1 .

RIPA lysis buffer containing protease inhibitor was used to extract tissue protein. Protein concentration was determined with BCA Protein Assay Kit. ColorMixed Protein Marker was applied as a protein size marker. Total proteins were separated by 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis and then transferred to a polyvinylidene difluoride (PVDF) membrane (3010040001, Roche, Germany). The PVDF membrane was then blocked with 5% milk for 2 h followed by incubation with the indicated primary antibody overnight at 4°C. After incubation with the corresponding secondary antibody at room temperature for 1.5 h, the protein bands were visualized using enhanced chemiluminescence (Tanon, China).

Tissues were fixed with 4% paraformaldehyde, dehydrated, embedded in paraffin and sectioned at 4μm. The sections were deparaffinized, then boiled with citrate solution, and incubated with 3% H₂O₂ for 10 min. After antigen retrieval and blocking, the slides were incubated with the indicated antibodies at 4°C overnight. Subsequently, the slides were incubated with the corresponding secondary antibody for 30 min at room temperature and then incubated with the streptavidin peroxidase complex. Staining was performed using a 3,3-diaminobenzidine (DAB) substrate kit for peroxidase reaction and counterstained with hematoxylin. Sections were washed, dehydrated, and sealed with neutral balsam. Finally, an Olympus microscope was used to acquire images.

All statistical analyses were conducted using Perl version 5.38.2 and R software version 4.3.1. A P-value of less than 0.05 was used to determine statistical significance. In MR analysis, we employed five methods, including “MR Egger,” “Weighted Median,” “Inverse Variance Weighting (IVW),” “simple mode,” and “weighted mode.” Among these, IVW is the primary method for causal inference due to its superior precision and stability. A P-value of less than 0.05 was considered indicative of a significant correlation between exposure and outcome. We conducted sensitivity and heterogeneity tests based on the Q statistic. The MR-egger regression test was utilized to assess the presence of horizontal pleiotropy based on its intercept term. MR-PRESSO compares the observed distance of all variants to the regression line with the expected distance under the null hypothesis of no horizontal pleiotropy. To evaluate whether the removal of individual influential SNPs influenced the overall estimate of causal effects, leave-one-out analyses were performed.

DEGs in the GEO database (GSE46602 and GSE79021) were identified in 189 cases of PCa and 63 normal controls using the R package “limma”. Due to batch effects present between different datasets (as shown in Figures 1A, C ), after removal, the distribution of data across different datasets became significantly more uniform, with medians aligned and both mean and variance showing stronger consistency(as shown in Figures 1B, D ). This indicates that after adjustment, the integrated data can be considered as a unified batch for subsequent analysis. Based on this, we identified a total of 21 DEGs in the training group, 11 of which were upregulated and 10 were downregulated. The expression levels of the 21 DEGs are displayed in heatmaps and volcano plots ( Figures 1E, F ).

We constructed a scale-free co-expression network using WGCNA to identify the most highly correlated modules in PCa. The optimal “soft” threshold β = 4 (scale-free R² = 0.9) was selected based on scale independence and mean connectivity ( Figures 2A, B ). A total of 5 modules with different colors were obtained through dynamic tree shearing ( Figure 2C ). Then, the modules were correlated with clinical characteristics to further identify the driver genes of the most strongly positively correlated module. We ultimately identified the ME turquoise module as being most highly correlated with DCM in PCa (R = 0.39, p < 0.001) ( Figure 2D ). A significant distribution of Gene Significance (GS) was observed across the five modules ( Figure 2E ). Therefore, we selected the turquoise module, which contains 1,195 genes, for further analysis. Additionally, we performed a correlation analysis between module membership (MM) and gene significance (GS), and found a significant positive correlation between them (R = 0.84, p < 0.001; Figure 2F ).

We intersected the DEGs with the WGCNA module genes to obtain a total of 20 genes ( Figure 3A ). To further assess the potential regulatory pathways of these genes in PCa, we conducted GO and KEGG enrichment analyses. The GO analysis indicated that these genes are primarily enriched in the biological process (BP) terms “intermediate filament organization,” “intermediate filament cytoskeleton organization,” and “intermediate filament-based process.” Regarding the cellular component (CC) ontology, these genes are involved in “intermediate filament” and “intermediate filament cytoskeleton.” The molecular function (MF) analysis showed that the genes are enriched in “structural constituent of cytoskeleton” and “scaffold protein binding” ( Figure 3B ). The KEGG enrichment analysis demonstrated that these genes participate in the “Staphylococcus aureus infection” and “Estrogen signaling pathway” ( Figure 3C ).

To construct a diagnostic model for PCa-related genes, we utilized a combination of 113 machine learning algorithms to analyze the previously selected 20 genes. In the training group dataset, we calculated that “RF” produced the best predictive performance among the 113 tested algorithms, and this model also performed well in the test group dataset. Furthermore, we calculated the AUC values for each model in both the training and validation groups, with the RF machine learning algorithm achieving the highest AUC value ( Figure 4A ). Therefore, we identified RF as the optimal diagnostic model. Additionally, we plotted the ROC curves for the RF model in the training and validation groups across three datasets ( Figures 4B–E ). The RF model identified six genes (SLC14A1, ARHGEF38, NEFH, MSMB, KRT23, KRT15) as the final validated core genes ( Figure 4F ).

To assess the diagnostic value of the six candidate core genes selected, we established ROC curves to evaluate the diagnostic specificity and sensitivity of each gene. The results are as follows: SLC14A1 (AUC 0.862), ARHGEF38 (AUC 0.817), NEFH (AUC 0.805), MSMB (AUC 0.784), KRT23 (AUC 0.728), KRT15 (AUC 0.824), all candidate genes have a high diagnostic value for PCa ( Figure 5A ). To delve deeper into the biological significance of these six core genes, we utilized GeneMANIA to analyze the relationships and interaction patterns between these genes and other closely related genes. we found 20 genes that are closely related to these 6 genes ( Figure 5B ). Next, we applied the Gene Set Variation Analysis (GSVA) method to analyze the signaling pathways in the KEGG database to assess the enrichment of each gene in these pathways. We found that the genes SLC14A1, MSMB, KRT23, and KRT15 are primarily enriched in “prostate cancer” signaling pathway ( Supplementary Figure S1 ). We also analyzed the expression level differences of these core genes in PCa patients and normal populations. The results showed that SLC14A1, NEFH, MSMB, KRT23, and KRT15 are underexpressed in PCa, while ARHGEF38 is highly expressed in PCa ( Figure 5C ). In addition, we also analyzed the correlations between these core genes and found that SLC14A1, NEFH, MSMB, KRT23, and KRT15 genes are all positively correlated with each other, while ARHGEF38 is negatively correlated with the other five genes ( Figure 5D ).

We used the CIBERSORT algorithm to evaluate the composition and abundance of 22 types of immune cell infiltration in the training group ( Figure 6A ). In addition, the correlation analysis of the 22 types of immune cells showed that Dendritic cells resting were positively correlated with Macrophages M1 (r = 0.62), T cells regulatory were positively correlated with NK cells activated (r = 0.50), while T cells CD4 memory resting were negatively correlated with T cells follicular helper (r = −0.70) ( Figure 6B ). The box plots showed that compared to the control group, PCa patients had higher levels of Plasma cells and Monocytes, while the levels of Mast cells resting were lower ( Figure 6C ).

To validate the relationship between the selected core genes and immune cell infiltration, we conducted a correlation analysis. The results indicated that SLC14A1 is positively correlated with Dendritic cells activated (R = 0.48, P = 0.027) and NK cells resting (R = 0.45, P = 0.039), and negatively correlated with Dendritic cells resting (R = -0.45, P = 0.039) and Macrophages M1 (R = -0.52, P = 0.016) ( Figure 7A ); ARHGEF38 is negatively correlated with T cells regulatory (Tregs) (R = -0.52, P = 0.017) and Mast cells resting (R = -0.56, P = 0.009) ( Figure 7B ); NEFH is positively correlated with NK cells activated (R = 0.46, P = 0.038) and negatively correlated with T cells CD4 memory activated (R = -0.58, P = 0.006) ( Figure 7C ); MSMB shows no significant correlation with immune cell infiltration ( Figure 7D ); KRT23 is positively correlated with Monocytes (R = 0.47, P = 0.031) and NK cells resting (R = 0.45, P = 0.040), and negatively correlated with Macrophages M1 (R = -0.45, P = 0.039), Dendritic cells resting (R = -0.48, P = 0.029), T cells CD4 memory activated (R = -0.49, P = 0.025), and B cells memory (R = -0.52, P = 0.015) ( Figure 7E ); KRT15 is positively correlated with NK cells resting (R = 0.48, P = 0.028) and negatively correlated with Macrophages M1 (R = -0.45, P = 0.039) ( Figure 7F ). Additionally, we visualized the correlations between the expression of all core genes and the infiltration of 22 immune cells ( Figure 7G ). These findings suggest that the associated core genes may contribute to PCa by influencing immune cell infiltration.

In our study, we initially conducted a two-sample MR analysis using eQTL data of six core genes combined with PCa data. Unfortunately, this analysis did not reveal any causal associations between these genes and PCa. Subsequently, we turned to another round of two-sample MR analysis using the pQTL data of these core genes along with PCa data. In this analysis, we found that only the MSMB pQTL had a significant causal relationship with the incidence of PCa (OR: 0.8589, 95% CI: 0.8085-0.9125, P < 0.05) ( Figure 8A ). This indicates that the expression of MSMB is a major protective factor in the pathogenesis of PCa, while the other core genes may only appear as biomarkers after the occurrence of PCa. The forest plot shows the effect sizes of various SNPs in the MSMB pQTL ( Figure 8B ). The scatter plot shows the regression trends of the five methods, all of which are consistent ( Figure 8C ). The funnel plot shows the distribution of causal effects ( Figure 8D ). The leave-one-out analysis reveals the impact of individual SNPs on the overall causal estimate ( Figure 8E ). The MR-Egger intercept analysis for the MSMB pQTL did not detect any potential horizontal pleiotropy (P = 0.55), indicating that the instrumental variables do not significantly affect the outcome through pathways other than the exposure. In the Cochran Q test for heterogeneity, heterogeneity was observed in the MSMB pQTL data (P = 0.013). This may be due to IVs originating from different analysis platforms, different studies or datasets, and different populations.

We measured the mRNA expression levels of six core genes in PCa tumor tissues and their corresponding para-cancerous tissues. Our results revealed that, compared to para-cancerous tissues, the expression levels of SLC14A1, NEFH, MSMB, KRT23, and KRT15 were significantly downregulated in PCa tumor tissues, while ARHGEF38 expression was notably upregulated, consistent with our initial predictions ( Figures 9A–F ). To further validate these findings at the protein level, we conducted Western blot analysis on tumor tissues and adjacent normal tissues from PCa patients. The results demonstrated that ARHGEF38 protein expression was significantly elevated in tumor tissues. In contrast, the expression levels of UT-B (encoded by SLC14A1), NEFH, MSMB, KRT23, and KRT15 were reduced in tumor tissues compared to adjacent normal tissues, although the decrease in UT-B did not reach statistical significance ( Figure 9G ). Quantitative analysis of the Western blot results is presented in Figure 9H , further supporting the observed protein expression trends.

Additionally, immunohistochemistry (IHC) was performed to evaluate the localization and expression levels of two key proteins in PCa and para-cancerous tissues. The IHC results revealed that ARHGEF38 protein was highly expressed in tumor tissues, while MSMB exhibited reduced expression in tumor tissues compared to para-cancerous tissues ( Figure 9I ). These findings were consistent with the mRNA and protein expression patterns observed in qRT-PCR and Western blot analyses.