Plasma proteome-based integrated mendelian randomization and multi-omics analysis identifies potential diagnostic biomarkers and therapeutic targets for prostate cancer

2.1 Data mining and collection

We obtained protein quantitative trait loci (pQTL) data from a study conducted by Ferkingstad et al. [15] using deCODE genetics (https://www.decode.com/summarydata/). The pQTL data satisfied the following criteria: (a) genome-wide significant associations (p < 5 × 10⁻⁸) and (b) independent associations (linkage disequilibrium (LD) clumping r² < 0.001). Summary statistics for PCa were acquired from the Integrative Epidemiology Unit (IEU) Open Genome-Wide Association Studies (GWAS) project (https://gwas.mrcieu.ac.uk/), with the analytical cohort (GWAS ID: ebi-a-GCST90018905, European) comprising 11,599 PCa cases and 199,628 controls, encompassing 211,227 single nucleotide polymorphisms (SNPs).

Besides, two bulk RNA-seq datasets (GSE200879 and GSE46602), a single-cell transcriptomic dataset (GSE193337), and a spatial transcriptomic dataset (GSE278936), containing human prostate cancer and normal prostate tissues, were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). The annotation information provided on the respective platforms was used as a reference for mapping corresponding gene symbols. Clinical information of patients with PCa and their corresponding mRNA expression data were downloaded from The Cancer Genome Atlas (TCGA) website (https://cancergenome.nih.gov/).

2.2 Mendelian randomization analysis

A bidirectional two-sample Mendelian randomization (MR) approach was employed to investigate the causal relationship between the human plasma proteome and prostate cancer (PCa). In this study, the plasma proteome was considered the exposure, and PCa was treated as the outcome, aiming to identify their potential roles as protective or risk-enhancing factors in PCa development.

Five analytical methods were applied to ensure robust causal inference, including inverse variance-weighted (IVW), MR-Egger, weighted median, simple mode, and weighted mode. IVW served as the primary method, while the other four methods were used as supplementary tools to validate the results. The IVW method combines Wald estimates from individual SNPs using a meta-analytic framework to generate an overall causal estimate, represented by the slope of the weighted regression of outcome effects on exposure effects [16]. The MR-Egger method was applied to assess horizontal pleiotropy, while the weighted median approach, recognized for its robustness in smaller datasets, served as a complementary method to enhance the reliability of the IVW and MR-Egger analyses [17]. The “TwoSampleMR” [18] and “MR-Pleiotropy RESidual Sum and Outlier (PRESSO)” [19] R packages were utilized in R (version 4.4.2) to perform the MR analyses and evaluate potential pleiotropy or other biases.

2.3 Sensitivity analysis and direction validation in MR

The MR-Egger intercept test [20] was applied to detect horizontal pleiotropy, with a statistically significant intercept indicating its presence. Cochran’s Q test [21] was used to evaluate SNP heterogeneity, and a significant Q statistic reflected the presence of heterogeneity. Sensitivity analysis was performed using the leave-one-out method, systematically excluding each SNP to assess the robustness of the results. Consistent MR estimates after the exclusion of individual SNPs demonstrated the reliability of the findings. A p-value < 0.05 was considered statistically significant. In the absence of pleiotropy and heterogeneity, the Inverse Variance-Weighted method was deemed valid.

2.4 Bulk RNA-seq data processing

Data normalization was performed using the limma package [22] to correct batch effects, with raw expression matrices log2-transformed and normalized via the NormalizeBetweenArrays function. Principal component analysis (PCA) confirmed the removal of batch effects. The limma package was used to identify differentially expressed genes (DEGs) with |logFC| >0.585 and an adjusted p-value < 0.05. Volcano plots and heatmaps of the top upregulated and downregulated genes were respectively generated using ggplot2 [23] and pheatmap packages [24], to visualize transcriptional changes.

2.5 Single-cell RNA-seq data processing

The single-cell RNA sequencing data were processed using standard workflows [25], including normalization, dimensionality reduction, clustering, cell type annotation, and intercellular communication analysis. A Seurat object was created, and low-quality cells (≤50 feature genes or ≥15% mitochondrial gene content) were removed. Data normalization was performed with the LogNormalize method (scale factor: 10,000), and highly variable genes were identified using variance-stabilizing transformation (VST). PCA was used for dimensionality reduction, with the top 20 components selected. Batch effects were corrected using the Harmony algorithm, and t-distributed stochastic neighbor embedding (t-SNE) was applied for visualizing cell distribution patterns.

Clustering analysis was conducted using the Louvain algorithm on shared nearest neighbor (SNN) graphs with a resolution of 0.6. Marker genes were identified through differential expression analysis (log2 FC >1, adjusted p-value < 0.05) and used to annotate cell clusters based on their expression patterns. Cell type annotation was further refined with the SingleR [26] tool and reference databases (e.g., HumanPrimaryCellAtlas and BlueprintEncode), identifying major cell types including cluster of differentiation (CD)4+ T cells, CD8+ T cells, B cells, monocytes, macrophages, fibroblasts, and endothelial cells. Manual curation and marker gene validation were conducted to ensure accurate annotation.

2.6 Spatial transcriptomic data processing

Spatial transcriptomic (ST) data were processed and analyzed to explore spatially resolved gene expression patterns [27]. The raw data were normalized and corrected using the SCTransform method, followed by dimensionality reduction (PCA), clustering, and visualization using Uniform Manifold Approximation and Projection (UMAP). The FindAllMarkers and FindMarkers functions, implemented in the Seurat R package, were used to identify DEGs derived from spatial transcriptomic data (ST-DEGs), which provide insights into spatially localized gene expression patterns and the molecular heterogeneity across tissue regions.

2.7 Multiple intercellular analysis

To investigate the spatial organization of cell types and their interactions within the tissue, we employed Multiple Intercellular Analysis (MIA), which integrates scRNA-seq data with ST data. Although ST data provide spatial resolution, they lack the cellular resolution of scRNA-seq. MIA overcomes this limitation by mapping scRNA-seq-derived cell type information onto the spatial context of ST data through shared gene expression profiles. This approach enables the identification of spatially resolved cell types, their functional states, and potential intercellular communication networks, providing a comprehensive view of the tissue microenvironment.

2.8 Functional enrichment analysis and PPI network construction

Gene Ontology (GO) analysis is a widely used method for annotating genes and their products by categorizing them into three main domains, including biological processes (BP), cellular components (CC), and molecular functions (MF). The Kyoto Encyclopedia of Genes and Genomes (KEGG) serves as a comprehensive resource that systematically analyzes biological pathways, integrating genomic, chemical, and functional information to elucidate molecular interactions and cellular processes. In this study, we utilized the DOSE (v3.22.0) [28] and clusterProfiler [29] (v4.4.4) packages in R (v4.2.1) to explore the biological relevance of the gene set and conduct pathway enrichment analyses. Furthermore, the enrichplot and ggplot2 packages were employed to generate high-quality visualizations of the results. The protein-protein interaction (PPI) network was constructed using the GeneMANIA [30] prediction server (http://www.genemania.org/).

2.9 Immune profiling analysis

The association between gene expression and the tumor immune microenvironment was analyzed using bulk RNA-seq data obtained from the TCGA-prostate adenocarcinoma (PRAD) cohort. Spearman correlation analysis was employed to assess the relationship between gene expression and immune cell infiltration levels, as well as the immune infiltration matrix data. Samples were stratified into high- and low-expression groups based on the median expression level of key biomarkers. Immune infiltration levels were quantified using the single-sample gene set enrichment analysis (ssGSEA) algorithm implemented in the Gene Set Variation Analysis (GSVA) package [31] (v1.46.0), with immune cell markers for 24 distinct immune cell types. The ggplot2 package was used to visualize the results.

2.10 Small molecule screening and molecular docking analysis

The relationships between small molecules and target proteins were obtained from the DsigDB [32] database (https://dsigdb.tanlab.org/). Network Venn diagrams were visualized using R software (v4.2.1) with the igraph (v1.4.1), ggraph (v2.1.0), and ggplot2 (v3.4.4) packages. Protein structure data were retrieved from the PDB database [33] (https://www.rcsb.org/), and the three-dimensional structures of small molecules were downloaded from the PubChem database [34] (https://pubchem.ncbi.nlm.nih.gov/). Molecular docking was performed using the CB-dock2 [35] online tool (https://cadd.labshare.cn/cb-dock2/index.php) to predict the binding modes between small molecules and target proteins.

2.11 Statistical analysis

Continuous variables were expressed as mean ± standard deviation (SD). For comparisons between two independent groups, the Student’s t-test was applied for data following a normal distribution, while the Wilcoxon rank-sum test was used for non-normally distributed data. Paired samples were analyzed using either the paired t-test or the Wilcoxon signed-rank test, based on the distribution of the data. Correlations between variables were assessed using the Pearson correlation coefficient. Diagnostic performance was evaluated using receiver operating characteristic (ROC) curves, and the area under the curve (AUC) was calculated to quantify accuracy. A significance threshold of p < 0.05 was applied for all statistical tests.

3.1 Flow chart and related data information

The workflow of this study design is illustrated in Fig. 1. The characteristics of each enrolled PCa dataset downloaded from GEO database are summarized in Table 1, and the detailed information is provided in Supplementary Table 1.

Fig. 1.The workflow of this study design. DEGs, differentially expressed genes; GO, Gene Ontology; GOLM1, Golgi membrane protein 1; IGFBP6, insulin-like growth factor binding protein 6; PCa, Prostate cancer; PRAD, prostate adenocarcinoma; ROC, receiver operating characteristic; SERPINF1, serpin family F member 1; ST, spatial transcriptomic; TCGA, The Cancer Genome Atlas; TGFB3, transforming growth factor beta 3; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.2 Identification of PCa-associated plasma proteins through MR analysis

After removing SNPs with insufficient instrument strength (F-statistic < 10), a total of 31,315 SNPs associated with 3457 plasma proteins were identified (Supplementary Table 2). To ensure the robustness of the instrumental variables (IVs), the distribution of F-statistics was assessed, ranging from 29.71 to 57,932.93, with a median of 54.75 and a mean of 184.27, demonstrating strong instrument strength across all included SNPs (Supplementary Fig. 1). MR analysis using these IVs identified 66 plasma proteins with IVW p-values < 0.05 (Supplementary Table 3). The circular plot shows the p-value distributions of plasma proteins across five MR methods (Fig. 2a). The volcano plot (Fig. 2b) highlights the regression coefficients (Beta values) and −log10(p-values) for these PCa-associated plasma proteins (PCa-PP), with significant proteins such as Collagen Type III Alpha 1 Chain (COL3A1) and Interleukin 1 Receptor Like 2 (IL1RL2) exhibiting significant associations. The Manhattan plot (Fig. 2c) further illustrates the genomic distribution of significant proteins across chromosomes, emphasizing the importance of COL3A1 and IL1RL2 in the causal relationship with PCa.

Fig. 2.Results of MR Analysis. (a) Circular plot showing the p-value distribution of plasma proteins across five MR methods. Each row represents a protein, and each column represents a method. (b) Volcano plot illustrating the distribution of Beta coefficients and −log10(p-values) for plasma proteins. (c) Manhattan plot displaying the genomic distribution of significant plasma proteins. Beta, Regression Coefficient; IVW Estimation, Inverse Variance Weighted Estimation; MR, Mendelian Randomization; Chr, Chromosome.

3.3 Identification of DEGs based on bulk RNA-seq data

After normalizing gene expression across datasets (Fig. 3a,b), 966 DEGs derived from bulk RNA-seq data (bulk-DEGs) were identified, including 335 upregulated and 631 downregulated genes (Fig. 3c, Supplementary Table 4). The top 20 upregulated and downregulated genes ranked by fold change were displayed in a heatmap, which showed significant differences in gene expression patterns between PCa and normal prostate tissue samples (Fig. 3d).

Fig. 3.Identification of DEGs based on bulk RNA-seq data. (a) Batch effect correction across datasets. (b) Principal component analysis (PCA) before and after batch effect correction. (c) Volcano plot of DEGs. Upregulated genes are shown in red, downregulated genes in green, and non-significant genes in gray. (d) Heatmap of the top 20 upregulated and downregulated DEGs. Red and blue indicate high and low expression levels by logFC, respectively.

3.4 Integration of single-cell RNA-seq and spatial transcriptomics

Single-cell RNA sequencing (scRNA-seq) identified 25 clusters, representing 11 major cell types, including immune cells (CD4+ T cells, CD8+ T cells, macrophages), fibroblasts, and epithelial cells (Fig. 4a). Differences in cell distribution were observed between normal and tumor samples. In tumor samples, CD4+ and CD8+ T cells exhibited a more localized and concentrated distribution compared to normal tissues, whereas epithelial and fibroblast cells showed a more dispersed distribution (Fig. 4b). Spatial transcriptomics (ST) analysis further identified 24 spatially resolved clusters, underscoring the cellular heterogeneity across samples (Fig. 4c). The DEGs associated with spatial transcriptomics data (ST-DEGs) between normal and tumor samples were identified for further analysis (Supplementary Table 5). Annotation of major clusters using MIA revealed that clusters 6 and 24 were significantly enriched with immune cells, including CD4+ T cells, CD8+ T cells, macrophages, and monocytes. In contrast, other clusters showed limited immune cell types, suggesting a distinct immune cell aggregation pattern in these clusters (Fig. 4d, Supplementary Table 6). Additionally, spatial visualization of selected normal and tumor sample slices demonstrated a marked increase in the absolute abundance of clusters 6, 24, and 23 in tumor samples, whereas clusters 7, 9, and 11 exhibited a relatively higher proportion in normal samples (Fig. 4e). This pattern indicates a significant enrichment of immune cell-related clusters in the tumor microenvironment, while stromal cell-related clusters are more abundant in normal tissue.

Fig. 4.Integration of single-cell RNA-seq and spatial transcriptomics. (a) Identification of 25 distinct cell clusters and 11 major cell types from the scRNA-seq dataset GSE193337. (b) Comparison of cell distribution patterns between normal and tumor samples. (c) Spatial transcriptomics analysis identified 25 spatially resolved clusters from the ST dataset GSE278936, revealing cellular heterogeneity across samples. (d) Annotation of major clusters using MIA analysis. (e) Spatial visualization of normal and tumor sample slices.

3.5 GO and KEGG analysis of PCa-PP, Bulk-DEGs, and ST-DEGs

The GO and KEGG enrichment analyses provided insights into the biological processes and pathways associated with PCa across different datasets. For PCa-PP, significant enrichment was observed in terms related to extracellular matrix dynamics, including extracellular matrix structural constituents, collagen-containing extracellular matrix, and processes like wound healing and regulation of inflammatory responses (Fig. 5a, Supplementary Table 7). These results highlight the role of extracellular remodeling and immune regulation in PCa progression. For Bulk-DEGs, pathways and biological functions such as glutathione metabolism, the phosphatidylinositol 3-kinase-protein kinase B (PI3K-Akt) signaling pathway, glycosaminoglycan binding, and epithelial cell proliferation were significantly enriched (Fig. 5b, Supplementary Table 8), reflecting alterations in metabolic processes, cell signaling, and extracellular matrix interactions in PCa. Finally, ST-DEGs revealed enrichment in functions such as complement and coagulation cascades, cell adhesion molecules, and extracellular matrix organization (Fig. 5c, Supplementary Table 9), emphasizing the importance of immune modulation, cell adhesion, and structural reorganization in the tumor microenvironment. These findings collectively underscore the critical roles of extracellular matrix remodeling, immune regulation, and metabolic shifts in PCa pathogenesis.

Fig. 5.Results of GO and KEGG analysis. (a) GO and KEGG enrichment analysis of PCa-associated plasma proteins (PCa-PP), (b) DEGs derived from bulk RNA-seq data (bulk-DEGs), and (c) DEGs associated with spatial transcriptomics data (ST-DEGs). BP, biological processes; CC, cellular components; MF, molecular functions; KEGG, Kyoto Encyclopedia of Genes and Genomes; PI3K, Phosphatidylinositol 3-Kinase; Akt, Protein Kinase B.

3.6 Key biomarkers identification and enrichment analysis

To identify key biomarkers associated with PCa progression, we integrated the results of PCa-PP and DEGs, identifying four key biomarkers, including TGFB3, IGFBP6, GOLM1, and SERPINF1 (Fig. 6a). A PPI network was constructed to visualize the interactions among these biomarkers and related proteins (Fig. 6b). Subsequently, GO and KEGG enrichment analyses were performed, and tumor progression-related terms were displayed in lollipop plots (Fig. 6c, Supplementary Table 10). The results showed involvement in epithelial to mesenchymal transition, SMAD signaling regulation, cytokine binding, extracellular matrix components, and pathways such as growth factor binding and the Wnt signaling pathway.

Fig. 6.Key biomarkers identification and enrichment analysis. (a) Venn diagram integrating PCa-PP and DEGs to identify key biomarkers. (b) PPI network visualizing the interactions among the four key biomarkers and related proteins. (c) Results of GO and KEGG enrichment analyses of the key biomarkers, with tumor progression-related terms displayed. MAPK, mitogen-activated protein kinase; PCa-PP, prostate cancer-associated plasma proteins; ST, spatial transcriptomics; SMAD, mothers against decapentaplegic homolog; DEGs, differentially expressed genes.

3.7 Key biomarkers expression validation and diagnostic performance in PCa

The TCGA-PRAD cohort was analyzed to validate the expression patterns of the four key biomarkers in PCa. The results revealed that IGFBP6, TGFB3, and SERPINF1 were significantly downregulated in tumor tissues compared to normal tissues, whereas GOLM1 was upregulated (Fig. 7a). Paired sample analysis further confirmed these expression patterns with significant differences (Fig. 7b). ROC curve analysis demonstrated that GOLM1 exhibited the highest diagnostic performance (AUC = 0.859), followed by SERPINF1 (AUC = 0.762), TGFB3 (AUC = 0.728), and IGFBP6 (AUC = 0.706) (Fig. 7c). Furthermore, a diagnostic model combining the four factors achieved an AUC of 0.911 (95% CI: 0.869–0.953), indicating excellent diagnostic accuracy (Fig. 7d). External validation was performed using independent datasets GSE46602 and GSE200879 to assess the robustness and generalizability of the diagnostic model (Fig. 7e–h). The model achieved AUCs of 0.994 in GSE46602 and 0.974 in GSE200879, strongly supporting its high diagnostic performance and potential clinical utility.

Fig. 7.Key biomarkers expression validation and diagnostic performance in PCa. (a) Expression levels of key biomarkers in PCa and normal prostate tissues. (b) Paired sample analysis validating the differential expression of key biomarkers. Statistical significance is denoted by ** and *** for p < 0.01 and p < 0.001, respectively. (c) ROC curve analysis illustrating the diagnostic performance of individual key biomarker in TCGA-PRAD cohort. (d) A combined diagnostic model integrating the four key biomarkers in TCGA-PRAD cohort. (e) ROC curve analysis illustrating the diagnostic performance of individual key biomarker in the GSE46602. (f) A combined diagnostic model integrating the four key biomarkers in the GSE46602. (g) ROC curve analysis illustrating the diagnostic performance of individual key biomarker in the GSE200879. (h) A combined diagnostic model integrating the four key biomarkers in the GSE200879. TCGA, The Cancer Genome Atlas; AUC, area under the curve; TMP, Tumor Microenvironment Proteins; TRP, Tumor-Related Proteins; FRP, Fibrosis-Related Proteins; PRAD, Prostate Adenocarcinoma; IGFBP6, Insulin-Like Growth Factor Binding Protein 6; TGFB3, Transforming Growth Factor Beta 3; GOLM1, Golgi Membrane Protein 1; SERPINF1, Serpin Family F Member 1.

3.8 Spatial expression patterns of key biomarkers in PCa

Spatial transcriptomics data were analyzed to further evaluate the expression patterns of the key biomarkers, revealing distinct regional differences between normal prostate and PCa tissues. The results demonstrated that TGFB3, IGFBP6, and SERPINF1 were significantly downregulated in PCa tissues compared to normal tissues, while GOLM1 was notably upregulated (Fig. 8). This spatial visualization highlights the contrasting expression profiles of these factors and their potential roles in PCa progression.

Fig. 8.Spatial expression patterns of key biomarkers in PCa. The spatial transcriptomics data were analyzed to evaluate the expression patterns of (a) TGFB3, (b) IGFBP6, (c) GOLM1, (d) and SERPINF1 in normal prostate and PCa tissues. TGFB3, Transforming Growth Factor Beta 3; IGFBP6, Insulin-Like Growth Factor Binding Protein 6; GOLM1, Golgi Membrane Protein 1; SERPINF1, Serpin Family F Member 1.

3.9 Immune microenvironment characteristics of key biomarkers in PCa

Immune cell enrichment score analysis was performed to compare high and low expression groups of TGFB3, IGFBP6, GOLM1, and SERPINF1 (Fig. 9a–d), with their correlations to various immune cell types summarized in a heatmap (Fig. 9e). As shown, TGFB3 and IGFBP6 demonstrated positive correlations with most immune cell types, including T cells, B cells, macrophages, and dendritic cells, highlighting their roles in promoting immune cell infiltration. Similarly, SERPINF1 exhibited positive correlations with a broad range of immune cell types. In contrast, GOLM1 displayed significant negative correlations with specific immune cells, such as neutrophils, mast cells, and macrophages, suggesting its potential role in suppressing their infiltration. Furthermore, correlation analysis with immune checkpoint molecules revealed that TGFB3, IGFBP6, and SERPINF1 were strongly and positively correlated with checkpoint molecules such as PDCD1, CD274, and CTLA4, whereas GOLM1 showed significant negative correlations with PDCD1, CTLA4, LAG3, and HAVCR2 (Fig. 9f).

Fig. 9.Immune microenvironment characteristics of key biomarkers in PCa. Boxplots comparing immune cell enrichment scores between high and low expression groups of (a) TGFB3, (b) IGFBP6, (c) GOLM1, and (d) SERPINF1. (e) Heatmap illustrating the correlation between key biomarkers and immune cell types in the TCGA-PRAD cohort. Positive correlations are shown in red, while negative correlations are shown in blue. (f) Heatmap showing the correlation between key biomarkers and immune checkpoint genes in the TCGA-PRAD cohort. Statistical significance is denoted by *, **, and *** for p < 0.05, p < 0.01, and p < 0.001, respectively.

3.10 Network construction and molecular docking with small molecules

To explore potential small molecules capable of targeting key biomarkers, an interaction network was constructed between candidate compounds and target proteins (Fig. 10a). The analysis revealed that several small molecules could simultaneously target two core proteins, such as diethylstilbestrol, irinotecan, 7646-79-9, camptothecin, selenium, and vitamin E. Among these, diethylstilbestrol was identified as interacting with both TGFB3 and IGFBP6 and was selected for molecular docking analysis. The three-dimensional (3D) structures of human TGFB3 (PDB ID: 2PJY) and IGFBP6 (PDB ID: 1RMJ) were obtained for docking analysis. The results demonstrated that diethylstilbestrol binds specifically to the active sites of TGFB3 (Vina score = −5.5) and IGFBP6 (Vina score = −4.6), forming stable complexes (Fig. 10b). These findings highlight diethylstilbestrol as a promising candidate for modulating the activity of TGFB3 and IGFBP6, offering potential therapeutic implications.

Fig. 10.Network analysis and molecular docking with small molecules. (a) Interaction network of small molecules and key biomarkers. (b) Molecular docking results of selected small molecules with TGFB3 and IGFBP3. GOLM1, Golgi Membrane Protein 1; TGFB3, Transforming Growth Factor Beta 3; IGFBP6, Insulin-Like Growth Factor Binding Protein 6; SERPINF1, Serpin Family F Member 1.