Abstract Glioblastoma (GBM), one of the most aggressive brain tumors, has a 5-year survival rate of less than 5%. Current standard therapies, including surgery, radiotherapy, and temozolomide (TMZ) chemotherapy, are limited by drug resistance and the blood–brain barrier. Integrating expression quantitative trait loci (eQTL) and protein quantitative trait loci (pQTL) data has shown promise in uncovering disease mechanisms and therapeutic targets. This study combined eQTL and pQTL analyses to identify potential GBM-related genes and circulating plasma proteins for therapeutic exploration. Using transcriptomic data from The Cancer Genome Atlas (TCGA), we identified 2,528 differentially expressed genes, including GPX7 and CXCL10. eQTL-MR analysis identifies GBM-associated differentially expressed genes and constructs a protein–protein interaction (PPI) network.Integrating pQTL data from the deCODE database, pQTL-MR, and colocalization analyses validated the therapeutic potential of GPX7 and CXCL10.These findings provide new perspectives on GBM biology and suggest actionable targets for therapy. Despite limitations due to sample size and population-specific data, this study highlights GPX7 and CXCL10 as promising candidates for further investigation and lays the foundation for targeted GBM treatments. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-025-13979-3. Keywords: Glioblastoma, eQTL, pQTL, Mendelian randomization, GPX7, CXCL10, Drug target Introduction Glioblastoma (GBM) is a highly invasive brain tumor, constituting 15% to 20% of primary brain tumors in adults. It is associated with a poor prognosis, with a 5-year survival rate of less than 5% [[32]1]. GBM exhibits high heterogeneity, with significant differences in the genome, phenotype, and immune microenvironment among tumors in different patients. This makes single-treatment approaches often ineffective, necessitating more personalized therapeutic strategies [[33]2]. Temozolomide (TMZ), as an effective treatment for GBM, is also limited in its effectiveness due to resistance and the constraints of the blood–brain barrier. Additionally, the GBM microenvironment contains various immunosuppressive factors that inhibit effective immune cell infiltration, rendering the tumor more resistant to treatment. The complexity of this microenvironment makes it challenging for conventional therapies to completely eliminate the tumor [[34]3]. A key cellular subtype contributing to the characteristics of GBM is a rare population of cells known as glioblastoma stem cells (GSCs), which are self-renewing, highly tumorigenic stem cells [[35]4]. GSCs drive the malignant progression of GBM through multiple mechanisms, including sustained proliferation, invasion, stimulation of angiogenesis, suppression of immune responses, and the development of drug resistance [[36]5]. Furthermore, GBM is a highly immunosuppressive tumor, where immune evasion mechanisms undermine the efficacy of single-agent immunosuppressive therapies, making them insufficient to alleviate GBM effectively [[37]6]. Consequently, researchers are actively developing new drugs for GBM therapy, aiming to discover more effective medications and treatment strategies. Randomized controlled trials (RCTs) are considered the gold standard for causal inference in epidemiological research [[38]7]. However, due to limitations such as medical ethics, participant compliance, and poor reproducibility of results, conducting RCTs in practice is often challenging [[39]8, [40]9]. In contrast, observational studies are relatively simple and easier to conduct, but due to issues like reverse causality and potential confounding factors, this method has faced criticism in terms of GBM etiology inference [[41]10, [42]11]. In recent years, the application of Mendelian randomization (MR) has provided an effective approach to address these issues [[43]12, [44]13]. The concept of MR is based on Mendel’s laws of inheritance, where parental alleles are randomly assigned to offspring, analogous to the randomization process in RCTs. MR utilizes genetic variations that are strongly associated with an exposure factor as instrumental variables (IVs) to infer the causal effect between the exposure and the outcome [[45]14]. Since genetic variations are innate, unaffected by environmental factors, and their association with the outcome follows a temporal causal sequence, MR analysis has become a prominent method for causal inference [[46]15]. In recent years, with the rapid advancement of genome-wide association studies (GWAS) and various omics data, the public sharing of large-scale summary statistics has provided new opportunities for the widespread application of MR in exploring the causal relationships between complex exposure factors and disease outcomes [[47]16]. In recent years, research combining expression quantitative trait loci (eQTL) and protein quantitative trait loci (pQTL) data has significantly established causal relationships between genes and diseases. This approach leverages genetic variations in gene expression and protein levels to provide deeper insights into how gene and protein functions jointly contribute to disease onset and progression. For example, these studies have proposed eQTL or pQTL as functional intermediates for investigating the potential biological mechanisms of neurodegenerative diseases [[48]17–[49]19]. Additionally, the joint analysis of eQTL and pQTL has provided more accurate guidance for identifying disease-related proteins. A recent study integrated pQTL and eQTL data, identifying PNKP as a therapeutic target and offering mechanistic insights into migraine pathophysiology [[50]20]. Proteins play a crucial role in various biological processes and are a major class of drug targets [[51]21]. Research has shown that disease-related protein drug targets supported by genetic associations are twice as likely to receive market approval compared to other protein drug targets [[52]22]. Plasma-circulating proteins have garnered increasing attention for their potential in GBM treatment, especially in diagnostics and therapeutic strategies. These proteins can serve as non-invasive biomarkers, aiding in early detection and monitoring of tumor progression [[53]23]. Circulating proteins reflect metabolic and immunological changes in GBM, which is critical for personalized treatment. For instance, proteins related to lipid metabolism, such as fatty acid synthase (FASN) and sterol regulatory element-binding protein (SREBP), have been linked to chemotherapy resistance. Targeting these proteins may enhance sensitivity to treatments such as TMZ [[54]24, [55]25]. Therefore, identifying and discovering plasma circulating proteins that are causally associated with GBM pathogenesis can provide real-time tumor biological information without the need for invasive tissue biopsies. This approach has the potential to significantly improve GBM management, aiding in early diagnosis, monitoring therapeutic response, and potentially identifying new therapeutic targets. MR analysis has become widely used in drug target development and drug repurposing [[56]26]. MR is a genetic instrumental variable analysis that typically utilizes single nucleotide polymorphisms (SNPs) from GWAS as genetic instruments to estimate the causal effect of exposures on outcomes. Compared to observational studies, MR can avoid the influence of confounding factors. With the advancement of high-throughput genomics and proteomics technologies in plasma, MR-based strategies have facilitated the identification of potential therapeutic targets for many diseases, such as stroke and Alzheimer's disease [[57]27, [58]28]. Although previous studies have successfully applied the combined MR analysis of eQTL and pQTL to identify cancer treatment targets, comprehensive analyses of this kind are still lacking in GBM research. Therefore, this study attempts to employ this innovative approach to identify potential therapeutic targets associated with GBM. In summary, this study aims to identify pathogenic genes causally associated with GBM by integrating eQTL and pQTL data and exploring the potential therapeutic value of plasma circulating proteins in GBM. We first utilized transcriptomic data from GBM and control groups in the TCGA database to screen for differentially expressed genes and performed eQTL-MR analysis to identify pathogenic genes associated with GBM. By constructing a protein interaction network and integrating plasma protein pQTL data from the deCODE database, we further validated and identified potential plasma protein targets for GBM through pQTL-MR analysis and colocalization analysis. This approach provides a new perspective for GBM treatment and may lay the foundation for the development of future targeted therapies. Materials & methods Study design The overall design of our study is depicted in Fig. [59]1. Initially, we analyzed transcriptomic data from GBM and control groups in the TCGA database to identify differentially expressed genes (DEGs). These DEGs were treated as exposures, with GBM serving as the outcome, to identify loci with potential causal associations to GBM. To validate the directionality of these associations, we conducted a bidirectional MRanalysis by reversing the roles of GBM and the identified loci, using GBM as the exposure and the loci as outcomes, to exclude the possibility of bidirectional causation. The MR analysis adhered to three core assumptions: (i) the Correlation Assumption, which requires a strong association between DEGs and instrumental variables; (ii) the Independence Assumption, ensuring no confounding factors are correlated with both the instrumental variables and GBM; and (iii) the Exclusivity Assumption, stipulating that the instrumental variables influence outcomes only through the exposure. Subsequently, we constructed a protein–protein interaction (PPI) network and integrated the plasma protein pQTL data from the deCODE database. We then performed pQTL-MR analysis to identify potential plasma proteins related to GBM as drug targets. To ensure the robustness of causal inference, we employed a range of MR methods. First, the F-statistic was calculated to evaluate the correlation between SNPs and the exposure, to screen for strong instrumental variables. Subsequently, the causal effect was estimated using the inverse variance weighting (IVW) method and Wald ratio estimation, while horizontal pleiotropy was assessed via MR-EGGER regression and the weighted median method. Heterogeneity was tested using Cochran’s Q statistic, and a leave-one-out analysis was performed to assess the impact of individual SNPs on the overall estimate. Additionally, Bayesian colocalization analysis was conducted to verify whether the observed associations are driven by the same causal variant rather than by linkage disequilibrium (LD) effects. Fig. 1. [60]Fig. 1 [61]Open in a new tab The overall flow chart of this study Data source and preparation The GBM dataset was procured from TCGA ([62]https://portal.gdc.cancer.gov/). This dataset comprised 170 GBM samples and 5 normal samples from the control group, from which their transcriptomic data were extracted. The datasets for eQTLs and GBM were procured from the Integrative Epidemiology Unit (IEU) database ([63]https://gwas.mrcieu.ac.uk/). GBM dataset originates from a European population and includes 91 cases and 174,006 controls, with a total of 16,380,303 SNPs. The data for patient prognosis analysis was obtained from the CGGA database. The dataset for pQTL-MR analysis in GBM was procured from the deCODE database. This dataset encompasses data from a population of 35,559 individuals in Iceland, with genome-wide association analysis conducted for 4,907 plasma circulating proteins. Analysis of differential gene in GBM The raw RNA-seq data from the TCGA-GBM dataset were processed using the R package DESeq2 ([64]https://bioconductor.org/packages/release/bioc/html/DESeq2.html) to identify DEGs. Genes with an adjusted p value < 0.05 and an absolute |log2FC|> 3 were deemed significantly differentially expressed. A volcano plot was then created to display the DEGs. Gene ontology, KEGG pathway, and GSEA analysis The obtained DEGs were subjected to enrichment analysis for Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and Gene Set Enrichment Analysis (GSEA) using the “ClusterProfiler” package. Enrichment with a p value < 0.05 was considered statistically significant. eQTL-MR analysis Bonferroni correction was applied to adjust for multiple testing, and eQTL SNPs with p < 5 × 10⁻⁸ and F > 10 thresholds were included in the MR analysis. DEGs identified through differential analysis were used as exposures, and GBM was set as the outcome. MR analysis was performed using the “TwoSampleMR” package ([65]https://github.com/MRCIEU/TwoSampleMR). For genes with only one available eQTL, the Wald ratio method was used, whereas for those with two or more genetic instruments, inverse-variance weighted (IVW) was applied to assess the robustness of the results. Odds ratio (OR) was used to reflect the relationship between gene expression and GBM risk: OR > 1 indicated that the gene was associated with an increased risk of GBM, OR < 1 suggested a potential protective effect of the gene, and OR = 1 implied no significant association. Utilize data from the CGGA database to retrieve and analyze the identified positive loci, clarifying the relationship between these genes and the prognosis of glioma patients. Protein-protein interaction network analysis After identifying positive loci with causal associations to GBM risk through the above steps, PPI analysis was performed for these loci. All PPI analyses were conducted using the GeneMANIA tool ([66]https://genemania.org/), and the results were visualized using Cytoscape. Detection of bidirectional causal relationship Applying the same eQTL selection criteria, bidirectional MR analysis was conducted on the dataset provided by the FinnGen database ([67]https://www.finngen.fi/fi) to detect potential reverse causal relationships. The MR-IVW method was utilized to estimate effects, ensuring the directionality of the association between positive genes and GBM prognosis. Results with a p value < 0.05 were considered statistically significant. pQTL-MR analysis To further explore the relationships between genes, proteins, and GBM, we extend the eQTL analysis by incorporating plasma cis-pQTL analysis. By examining the association between plasma protein expression levels and GBM risk, we aim to identify gene regulatory mechanisms that might drive protein level changes. The 29 positive loci associated with increased GBM risk were queried in the deCODE dataset, and the results were analyzed using the “TwoSampleMR” package. For proteins with only one available pQTL, the Wald ratio method was used, while for those with two or more genetic instruments, IVW was applied to assess the robustness of the findings. OR were used to reflect the relationship between gene expression and GBM risk: OR > 1 indicated that the gene was associated with an increased risk of GBM, OR < 1 suggested a potential protective effect of the gene, and OR = 1 implied no significant association. Sensitivity analysis and bayesian co-localization analysis Heterogeneity was evaluated using Cochran's Q test, and when heterogeneity was present, the IVW random-effects model was utilized [[68]29]. The Leave-One-Out approach assessed the influence of excluding each SNP on the outcome. The MR-EGGER method, with an intercept test, was used to detect horizontal pleiotropy. To investigate whether two correlated signals (protein and GBM risk) share a common causal variant while distinguishing linkage disequilibrium (LD) confounding, Bayesian co-localization analysis was conducted based on summary statistics from GWAS of protein and GBM, using the "coloc" package [[69]30]. The analysis considered five hypotheses: (i) no causal variant exists in the genomic locus (H0) for either protein or GBM; (ii) only protein (H1) has a causal variant; (iii) only GBM (H2) has a causal variant; (iv) protein and GBM (H3) have distinct causal variants; (v) protein and GBM (H4) share a common causal variant. For each protein, SNPs within ± 1Mb of the pQTL were included. When a protein had multiple pQTLs, co-localization analysis was conducted separately for each, focusing on the pQTL showing the strongest co-localization evidence. Default parameters were used, with p1 = 1 × 10^−4 (prior probability of an SNP being associated with the protein), p2 = 1 × 10^−4 (prior probability of an SNP being associated with GBM), and p12 = 1 × 10^−5 (prior probability of an SNP being associated with both protein and GBM). Given the sensitivity of co-localization to priors and window size, additional analyses were performed with alternative priors (p12 = 1e − 6) and a window size of ± 250 kb to assess robustness. Posterior probabilities were used to assess support for each hypothesis. A posterior probability greater than 80% for H4 (PP4) under different priors and window sizes was considered strong evidence for colocalization [[70]30–[71]32]. Visualization of co-localization results was conducted using the "LocusCompareR" package [[72]33]. All analyses were performed using R version 4.2.1. Results Selection of significantly DEGs in GBM RNA-Seq data from the TCGA database were processed using R, resulting in the identification of 23,557 DEGs. These genes were further filtered, with a threshold of adjusted P < 0.05 and |log2FC|> 3 set for significance. After filtering, 2,528 significant differentially expressed genes met the criteria and were included in the subsequent analysis. Among these, 1,440 genes, including GPX7 and CXCL10, were upregulated, while 1,088 genes, including TRIM17 and SLC7A4, were downregulated (Fig. [73]2A). Fig. 2. [74]Fig. 2 [75]Open in a new tab Differential gene screening and enrichment analysis of GBM and normal samples. A Volcano diagram of DEGs in GBM and normal samples. B GO enrichment analysis of DEGs. C KEGG enrichment analysis of DEGs. D GSEA enrichment analysis of DEGs GO, KEGG, and GSEA analysis of DEGs We conducted enrichment analysis, including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA), on the 2,528 DEGs. GO analysis, which included Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF), highlighted cellular-level changes. BP analysis revealed modifications in cross-synaptic signaling involved in chemical synaptic transmission, vesicle-mediated transport in the synapse, synaptic vesicle cycle, and neurotransmitter transport. CC analysis showed alterations in the synaptic membrane, neuron cell body, postsynaptic membrane, and ion channel complexes. MF analysis pointed to dysregulation in channel activity, passive transmembrane transporter activity, ion channel activity, gated channel activity, and voltage-gated ion channel activity. These results suggest that GBM alterations primarily occur at the cellular level, affecting synaptic regulation, neurotransmitter transport, and channel activity at both biological and molecular levels (Fig. [76]2B).KEGG pathway enrichment analysis indicated that differentially expressed genes are involved in pathways related to neuroactive ligand-receptor interactions, the calcium signaling pathway, the cAMP signaling pathway, morphine addiction, GABAergic synapse, cholinergic synapse, serotonergic synapse, synaptic vesicle cycle, taste transduction, and nicotine addiction (Fig. [77]2C).GSEA enrichment analysis highlighted significant enrichment in pathways associated with allograft rejection, autoimmune thyroid disease, cell cycle, graft-versus-host disease, ribosome, and systemic lupus erythematosus. These findings imply potential links between GBM and the identified pathways (Fig. [78]2D). Identification of positive loci causally associated with the onset of GBM The eQTLs derived from the above DEGs were used as exposures, with GBM as the outcome. MR analysis was performed using the “TwoSampleMR” package. The results revealed a total of 29 positive genes (p < 0.05). Among these, CLEC4E, LOC283194, AIDAP2, DDIAS, ACTN1-DT, SYT1, HNRNPA1P21, GBP1P1, IGLV3-9, SLC16A4, HLF, PHF24, MIR657, SLC7A4, IGLV2-14, ST13P19, and TRIM17 were associated with a decreased risk of GBM, while IGLV1-47, CDC45, CREB3L3, STX1A, PITPNM3, IGLV4-69, XRCC6P2, RAPGEFL1, CABP1, CXCL10, GPX7, and HOXB2 were associated with an increased risk of GBM (Fig. [79]3A, B and Supplementary Material 1). These above 29 positive loci were selected for further analysis. The data analysis from the CGGA database identified 15 available genes, most of which have expression levels associated with the prognosis of glioma patients (Supplementary Material 2). Fig. 3. [80]Fig. 3 [81]Open in a new tab eQTL-MR analysis and PPI network analysis. A The OR values of positive loci were identified from the MR analysis, where OR reflects the relationship between gene expression and GBM risk. An OR > 1 indicates a positive association with GBM, OR < 1 suggests a protective role, and OR = 1 indicates no significant association. B Results of the eQTL-MR analysis, where the horizontal axis represents the MR estimated causal effect values of positive loci on GBM, and the vertical axis lists the gene names. The color of each point represents the p-value from the MR analysis, and the size of each point corresponds to the absolute value of the MR estimate. C A PPI network was constructed from GENEMANIA using the 17 core genes identified through the analysis. The network illustrates the protein–protein interactions between these genes, highlighting their potential roles in GBM pathogenesis PPI network of positive loci causally associated with the onset of GBM A total of 29 positive loci were imported into GENEMANIA to construct a PPI network, identifying 17 key genes, including GPX7 and CXCL10. However, genes such as MIR657, IGLV4-69, IGLV1-47, IGLV2-14, IGLV3-9, GBP1P1, ST13P19, HNRNPA1P21, XRCC6P2, AIDAP2, LOC283194, and ACTN1-DT were not successfully identified. The generated results were exported to Cytoscape for PPI network visualization. The relationships among the key genes are shown in the figure. Based on the GeneMANIA tool, a gene network containing various types of interactions, including co-expression, physical interactions, predicted interactions, co-localization, shared protein domains, genetic interactions, and pathways, was constructed for target genes such as CXCL10 and GPX7. The analysis revealed that these genes and their associated genes were significantly enriched in multiple biological functions, including neurotransmitter release, synaptic vesicle-mediated signaling, neurotransmitter transport regulation, and synaptic cycling. The network weights were calculated using an automatically selected weighting method, with data sourced from multiple high-quality databases and research literature. Specifically, the analysis showed that co-expression networks accounted for the highest proportion (44.40%), including studies related to genes in the human central nervous system and those associated with liver expression quantitative trait loci. The physical interaction network accounted for 12.90%, with primary data derived from the Reactome pathway database and BioGRID studies. The co-localization network accounted for 11.55%, revealing the spatial relationships of genes within the cell, suggesting that they may function synergistically within the same functional modules. Furthermore, pathway analysis indicated that these genes are involved in several key biological processes, including immune responses, inflammatory regulation, oxidative stress balance, and extracellular matrix remodeling. CXCL10 was found to be a core node in immune regulation and inflammation pathways, closely related to immune cell recruitment and tumor microenvironment regulation. GPX7 plays a critical role in redox balance and DNA damage repair, forming a significant interaction network with antioxidant-related genes. These results suggest that the target genes may play an important role in tumor immune microenvironment regulation and neurofunctional control, providing data support for further exploration of their molecular mechanisms and potential therapeutic targets. (Fig. [82]3C and Supplementary Material 3). Bidirectional MR analysis The 29 positive loci obtained were subjected to bidirectional MR validation using data from the FinnGen database. GBM was considered as the exposure factor, and the eQTLs of DEGs were used as the outcome. The results showed that, except for HOXB2 (p < 0.05), the outcomes of bidirectional MR for the remaining 28 positive loci were not statistically significant (p > 0.05). This indicates that there is a unidirectional causal relationship between the 28 positive loci and GBM. There may be a bidirectional causal relationship between HOXB2 and GBM (Table [83]1). Table 1. Results of bidirectional mendelian randomization Outcome Exposure Method nsnp b se pval CREB3L3 Glioblastoma (age-stratified), Glioblastoma IVW 5 -0.007106561 0.015043951 0.636650852 7 -0.006452276 0.014307202 0.65200327 SYT1 5 -0.016673008 0.015044119 0.267743441 7 -0.006931251 0.014307641 0.628070769 PITPNM3 4 -0.00282299 0.02082459 0.892168805 5 0.018330425 0.016442027 0.264913895 CDC45 5 0.008613326 0.015044088 0.566956949 7 0.005363256 0.014307527 0.707767914 SLC7A4 5 0.004073113 0.015044236 0.786589114 7 3.48E-05 0.014307567 0.998060515 STX1A 5 0.002211464 0.018281167 0.903715189 7 -0.002847035 0.018229146 0.875890822 RAPGEFL1 5 -0.012778424 0.016518757 0.439184795 7 0.007885057 0.018702186 0.673308619 HLF 5 0.003099699 0.015044225 0.836760355 7 0.001515787 0.014307551 0.915627457 GPX7 5 0.005373937 0.015206469 0.72379001 7 0.001340018 0.014933143 0.928498207 PHF24 5 -0.000358246 0.01504424 0.9810019 7 0.005463229 0.01430756 0.702578504 CABP1 5 0.014991238 0.020739311 0.469777183 7 0.002943062 0.014307281 0.837022267 TRIM17 5 0.024295665 0.021592961 0.260518719 7 0.040909871 0.026536513 0.123159995 DDIAS 5 -0.018642862 0.015044052 0.215264719 8 -0.026084294 0.01411218 0.064551255 CLEC4E 5 -0.010331365 0.015044072 0.492246644 7 -0.000867175 0.01430737 0.951669539 SLC16A4 5 -0.00382559 0.019239336 0.842386244 7 -0.015755651 0.014307322 0.270796599 CXCL10 5 -0.002740918 0.019889262 0.890391422 7 -0.026586929 0.014306995 0.063124042 HOXB2 5 0.021781401 0.015043922 0.147657918 7 0.036168728 0.014544894 0.012893556 MIR657 5 -0.00551235 0.022586993 0.807192295 7 -0.009268314 0.016464526 0.573485094 IGLV4-69 5 0.01182686 0.01504401 0.431779176 7 0.007468506 0.01461845 0.609423958 IGLV1-47 5 -0.004983835 0.015044053 0.740430917 7 0.000164123 0.014307441 0.990847536 IGLV2- 14 5 0.005923679 0.015044251 0.693765354 7 0.01039328 0.014307517 0.467580725 IGLV3-9 5 -0.012165302 0.015044195 0.418723612 7 -0.01723212 0.014307504 0.228430622 GBP1P1 5 0.014555461 0.016705541 0.3835929 7 0.020306347 0.014307184 0.155808296 ST13P19 5 0.011086437 0.017395695 0.523923446 7 0.01013621 0.014307121 0.478651534 HNRNPA1P21 5 -0.018392867 0.016442688 0.263308824 7 -0.020117034 0.014307172 0.159700254 XRCC6P2 5 -0.027347765 0.015043419 0.06907602 7 -0.024392116 0.014306977 0.088211113 AIDAP2 5 0.014467611 0.015044038 0.336208352 7 0.012005673 0.014307392 0.401399874 LOC283194 5 0.008283151 0.015044063 0.581912952 7 0.007804019 0.014307461 0.585443369 ACTN1-DT 5 0.007656666 0.015044149 0.61078979 7 0.016712846 0.014307402 0.242756016 [84]Open in a new tab Identification of plasma circulating proteins causally associated with the onset of GBM To identify plasma circulating proteins potentially involved in the onset of GBM, we queried the 29 significant genetic loci identified in our eQTL-MR analysis against the deCODE database, which provides pQTLs associated with circulating protein levels. From this query, we identified three available plasma proteins: GPX7, CXCL10, and STX1A. For inclusion in our MR analysis, we applied a stringent selection criterion: only SNPs that meet the significance threshold (p < 5 × 10⁻⁸, F > 10) were included. Additionally, we required that multiple independent pQTLs be available for robust instrumental variable construction. Our analysis revealed that GPX7 and CXCL10 met these criteria and were significantly associated with an increased risk of GBM (Fig. [85]4, Tables [86]2 and [87]3). Specifically, GPX7 demonstrated strong associations with loci previously implicated in oxidative stress regulation, while CXCL10, a well-known chemokine involved in immune modulation, exhibited significant links to inflammatory pathways relevant to glioblastoma progression. In contrast, STX1A, despite being identified in our database query, had only one available SNP that did not meet our inclusion criteria for pQTL-based MR analysis. Consequently, STX1A was excluded from further causal inference analysis. To further validate these findings, we will conduct sensitivity analyses, including MR-Egger regression and weighted median estimation, to assess the robustness of causal associations and potential pleiotropy. Additionally, we will compare our results with existing literature on GPX7 and CXCL10 in glioblastoma to contextualize their potential biological relevance. Fig. 4. [88]Fig. 4 [89]Open in a new tab pQTL-MR analysis and sensitivity analysis of plasma circulating proteins associated with GBM onset. A MR analysis of GPX7. B The relationship between GPX7 and the SNP effect of GBM. With the increase of the SNP effect of GPX7, the SNP effect of GBM is also increasing. C MR analysis results of CXCL10. D The relationship between CXCL10 and the SNP effect of GBM. As the SNP effect of CXCL10 increases, the SNP effect of GBM also increases Table 2. Results from two sample MR of GPX7 Outcome Exposure Methods nSNP b se pval GBM GPX7 MR Egger 1325 0.3844 0.1132 0.0007 Weighted median 0.3077 0.1084 0.0045 Inverse variance weighted 0.1830 0.0633 0.0038 Simple mode -0.0430 0.3298 0.8965 Weighted mode 0.4933 0.2083 0.0180 [90]Open in a new tab Table 3. Results from two sample MR of CXCL10 Outcome Exposure Methods nSNP b se pval GBM CXCL10 MR Egger 387 4.8031 0.8337 1.7115E-08 Weighted median 3.8828 0.1768 6.9545E-107 Inverse variance weighted 3.3455 0.1275 8.4279E-152 Simple mode 3.8954 0.4934 3.0287E-14 Weighted mode 3.8954 0.4760 4.5306E-15 [91]Open in a new tab Sensitivity analysis and colocalization results for GPX7 and CXCL10 Sensitivity analysis for GPX7 and CXCL10 indicated no significant heterogeneity in the results (p > 0.05) (Supplementary Material 4). However, horizontal pleiotropy was detected for GPX7 (p < 0.05), which violates a core assumption that instrumental variables affect outcomes solely through exposure. This finding implies that certain SNPs associated with GPX7 might influence GBM risk not only through GPX7 expression but also via other unknown biological pathways. Such pleiotropic effects could confound the causal inference between GPX7 and GBM. To mitigate this problem, we employed MR-Egger regression to quantify the bias introduced by horizontal pleiotropy and to provide adjusted causal effect estimates [[92]34]. Despite the detected pleiotropy, MR-Egger yielded robust causal estimates, further supporting the causal link between GPX7 and GBM. A "Leave-one-out" sensitivity analysis was conducted, progressively removing individual SNPs and recalculating the meta-effect for the remaining SNPs. The results showed consistent overall error lines, which remained to the right of zero, demonstrating that horizontal pleiotropy had a limited impact on the findings (Supplementary Material 4). To further validate the relationship between GPX7 and CXCL10 with, GBM, colocalization analysis was performed to assess whether signals for protein levels and GBM risk originated from shared causal variants. Results revealed a shared causal variant between GPX7 and GBM, with a posterior probability exceeding 80%. Specifically, colocalization was identified at the rs141755176 locus for GPX7's pQTL and GBM GWAS (Fig. [93]5A). CXCL10 showed multiple colocalized loci with similarly strong evidence (Fig. [94]5B). This colocalization analysis effectively ruled out false-positive associations due to LD, bolstering the reliability of the causal inference. Fig. 5. [95]Fig. 5 [96]Open in a new tab The results of co-localization analysis of GPX7 and CXCL10 with GBM. A The GPX7-pQTL and GBM-GWAS were successfully co-located, and the two shared a common variant. B The CXCL10-pQTL and GBM-GWAS were successfully co-located, and the two shared a common variant Discussion In this study, we integrated eQTL and pQTL data to investigate pathogenic genes associated with GBM and the potential therapeutic value of circulating plasma proteins. Our findings indicate that GPX7 and CXCL10 exhibit significant causal relationships with GBM pathogenesis and may serve as potential therapeutic targets. Since the human proteome represents a major source of therapeutic targets for tumors, we employed a comprehensive analytical approach to identify potential drug targets for GBM. Initially, we screened DEGs between control tissue and GBM tissue by analyzing RNA-Seq data from TCGA, identifying 2,528 significant DEGs that may play crucial roles in GBM pathogenesis. To further elucidate the functions of these genes, we conducted extensive enrichment analyses, including GO, KEGG, and GSEA. GO analysis revealed that DEGs were predominantly enriched in biological processes related to synaptic signaling, neurotransmitter transport, and calcium signaling pathways, suggesting that GBM development may involve disruptions in neuronal signal transmission and synaptic dysfunction. Specifically, synaptic vesicle cycling, neurotransmitter transport, and trans-synaptic signaling were significantly affected, emphasizing the potential association between GBM onset and neuronal dysfunction. KEGG pathway analysis further demonstrated that these DEGs were involved in multiple neural function-related pathways, such as neuroactive ligand-receptor interactions, calcium signaling, cAMP signaling, and GABA receptor signaling, which may undergo substantial alterations in GBM. GSEA analysis identified significant enrichment in immune response, cell cycle regulation, and allograft rejection pathways, implying that GBM pathogenesis is not solely driven by tumor cell proliferation but also involves immune evasion and remodeling of the tumor microenvironment. In recent years, studies have suggested that neurotransmitters and synaptic activity may play a crucial role in GBM progression. For example, glutamate promotes GBM cell proliferation and invasion through AMPA and NMDA receptors and may reshape the tumor microenvironment via excitotoxicity [[97]35, [98]36]. Additionally, GABAergic signaling may have an inhibitory effect on tumor growth and contribute to tumor microenvironment regulation [[99]37]. Recent studies have also discovered that glioma cells may form functional synapses with neurons, thereby promoting tumor progression through neuro-glial interactions [[100]38]. To identify potential causal genes in GBM, we utilized eQTL data to screen genes with strong causal associations with GBM risk, identifying 29 positive loci significantly linked to GBM. CXCL10 plays a critical role in the tumor immune microenvironment, as it not only facilitates immune cell recruitment but also modulates immune responses, contributing to immune evasion. Upregulation of CXCL10 expression may promote immune escape, thereby exacerbating GBM progression. GPX7, on the other hand, is essential for oxidative stress regulation and DNA repair, particularly in tumor cells responding to oxidative damage and DNA repair mechanisms. Aberrant GPX7 expression may lead to redox imbalance, consequently promoting tumor cell growth and metastasis. The causal relationships identified through MR may be influenced by linkage disequilibrium (LD), reverse causation, horizontal pleiotropy, or genetic confounding [[101]21]. Therefore, we performed bidirectional MR analysis, which confirmed that, except for HOXB2, none of the identified genes exhibited evidence of reverse causality [[102]39]. To further explore interactions between GPX7, CXCL10, and other genes, we conducted a PPI network analysis, revealing that these genes closely interact with others involved in synaptic vesicle-mediated signaling and neurotransmitter transport. SYT1 overexpression has been shown to inhibit GBM cell proliferation and promote apoptosis [[103]40]. Additionally, SYT1 is mutated in GBM, suggesting its potential as a therapeutic target [[104]41]. PITPNM3 exhibits lower tumor-region expression compared to the peritumoral area in female GBM patients, indicating a possible sex-specific role [[105]42]. CDC45 regulates DNA replication through the CDC45-MCM-GINS (CMG) complex, influencing GBM cell proliferation and apoptosis [[106]43]. HLF suppresses GBM cell growth and enhances the efficacy of TMZ treatment [[107]44]. Increased expression of CABP1 is negatively correlated with progression-free survival in GBM patients [[108]45]. TRIM17 overexpression significantly inhibits GBM cell proliferation, while its silencing produces the opposite effect [[109]46]. Furthermore, SLC16A4 (MCT4) is highly expressed in GBM, particularly in the more invasive mesenchymal subtype. Its inhibition affects the HIF1α-related signaling pathway, suggesting its potential as a novel therapeutic target [[110]47]. Notably, STX1A is a key regulatory protein involved in synaptic vesicle membrane fusion. It primarily facilitates synaptic vesicle docking and neurotransmitter release by participating in SNARE complex formation [[111]48]. Additionally, STX1A may influence tumor cell energy metabolism by regulating the glucose transporter GLUT1, highlighting its potential role in GBM that warrants further investigation [[112]49]. For CREB3L3, RAPGEFL1, PHF24, CLEC4E, SLC7A4 and DDIAS, no direct studies have been found linking them to GBM. Their biological mechanisms in GBM require further exploration in future studies. Subsequently, we utilized MR analysis of published plasma circulating protein pQTL data to evaluate proteins with causal associations with GBM. We exclusively used cis-pQTLs of plasma proteins as instrumental variables to identify drug targets, as they directly influence transcription and/or translation [[113]50]. The results further confirmed that GPX7 and CXCL10 are potential therapeutic targets for GBM. Sensitivity analysis indicated that GPX7 exhibited horizontal pleiotropy, suggesting that certain SNPs within this gene might influence GBM risk through unknown biological pathways. We corrected for horizontal pleiotropy using MR-Egger regression, reinforcing the causal relationship between GPX7 and GBM. To mitigate biases arising from horizontal pleiotropy, Bayesian colocalization analysis was employed to exclude confounding effects introduced by LD. Colocalization analysis revealed that GPX7 and CXCL10 shared multiple causal variants with GBM, eliminating false positives due to LD and strengthening the robustness of our causal inference. However, these associations do not fully elucidate the relationship between the identified proteins and GBM. Thus, the findings should be interpreted cautiously. Human GPX7 is a nonselenocysteine-containing neutral antioxidant enzyme found in mammals, also known as nonselenocysteine phospholipid hydroperoxide GPX (phGPX = GPX4) due to its homology with phospholipid hydroperoxide GPX (GPX4). GPX7 serves as a crucial sensor of oxidative and endoplasmic reticulum (ER) stress [[114]51]. The redox state of GPX7 is determined by the cellular oxidative status. During oxidative stress, GPX7 becomes activated and transfers disulfide bonds to specific proteins, enhancing their activity to respond to oxidative challenges. Currently, limited research has explored the role of GPX7 in gliomas. The few available studies have reported that elevated GPX7 expression is associated with poor outcomes in gliomas. Silencing GPX7 has been shown to enhance ferroptosis-related oxidative stress in glioma cells [[115]52]. The upregulation of GPX7 is tightly regulated by epigenetic processes, which significantly affect the overall survival of patients with lower-grade gliomas (LGG) [[116]53]. Bioinformatics analyses indicate that GPX7 is upregulated across various tumor types and serves as a potential prognostic biomarker for LGG. Moreover, high GPX7 expression can predict the sensitivity of LGG patients to TMZ treatment [[117]54]. Studies have shown that the upregulation of GPX7 is tightly regulated by epigenetic mechanisms. In LGG patients, GPX7 is associated with processes such as immune responses and synaptic transmission, whereas in GBM, it is primarily linked to metabolic regulation. GPX7 is also closely related to immune cell infiltration and immune modulation within the tumor microenvironment, significantly influencing local immune responses [[118]53]. In our analysis, we identified confounding horizontal pleiotropy in GPX7. Horizontal pleiotropy refers to a genetic variant being associated with multiple phenotypes through different biological pathways [[119]55]. Specifically, SNPs may influence GBM progression without affecting GPX7 protein expression, primarily due to the diverse functions of these SNPs. For example, SNPs located in coding regions (exons) may lead to amino acid substitutions, directly impacting protein structure and function [[120]56]. SNPs may also affect DNA methylation, histone modification, and chromatin accessibility [[121]57]. Additionally, SNPs in the 3' untranslated region (UTR) could alter microRNA binding sites, thereby affecting mRNA degradation or translational regulation [[122]58]. Some SNPs are found in drug-metabolizing enzyme genes or drug target genes, influencing drug metabolism rates, efficacy, or toxicity (pharmacogenomics) [[123]59]. Therefore, caution is warranted when interpreting the relationship between GPX7 and GBM. Human CXCL10 is a 10 kDa protein functionally classified as a Th1-type chemokine. It binds to the CXCR3 receptor and regulates immune responses by activating and recruiting leukocytes such as T cells, eosinophils, and monocytes [[124]60]. In our study, CXCL10 emerged as another potential drug target. Current research on CXCL10 primarily focuses on tumor immunity. For instance, in colorectal cancer, CXCL10 expression has been associated with intratumoral CD8^+ T cell infiltration and the reprogramming of the tumor vascular system [[125]61]. In hepatocellular carcinoma, CXCL10 has been found to modulate the tumor microenvironment (TME) associated with fibrosis, thereby influencing tumor progression [[126]62]. In GBM-related studies, the therapeutic approach of upregulating CXCL10 expression in the TME has been recognized as a potential strategy. It may increase tumor-infiltrating T cells and enhance their activity. However, effective delivery methods are still lacking. Researchers have demonstrated that peritumoral administration of CXCL10 and Nrf2-overexpressing mesenchymal stem cells (MSCs) under MRI guidance can significantly restrict GBM growth by activating T lymphocytes within the TME [[127]63]. These findings underscore the potential of CXCL10 as a drug target. In our study, we observed that an increase in SNP effects associated with CXCL10 corresponded to an increase in SNP effects linked to GBM. However, research exploring the relationship between CXCL10 and GBM remains limited, marking this as a promising direction for future investigation. In this study, a bidirectional causal relationship was identified between HOXB2 (Homeobox B2) and GBM, suggesting that HOXB2 may act both as a contributing factor to GBM development and as a target influenced by GBM progression. HOXB2, a member of the HOX gene family, functions as a transcription factor critical for regulating cell differentiation, migration, and positioning during embryonic development. Aberrant expression of HOX genes has been implicated in various cancers, particularly in promoting cancer stem cell properties, invasiveness, and metastasis [[128]64, [129]65]. The overexpression of HOXB2 in GBM may enhance the proliferation and invasion of tumor cells, suggesting its role as a potential pathogenic factor in GBM [[130]66]. However, studies have shown that HOXB2 not only promotes tumor development but also plays a role in remodeling the extracellular matrix (ECM), thereby limiting the progression of triple-negative breast cancer [[131]67]. Therefore, the bidirectional causal relationship between HOXB2 and GBM may, in part, be attributed to its regulation of the ECM in GBM. The excessive ECM in GBM may, in turn, influence HOXB2 expression. In addition, the aberrant activation of numerous signaling pathways in GBM may contribute to the regulation of HOXB2 expression by GBM. For example, previous studies have demonstrated that the Wnt signaling pathway is abnormally activated in GBM [[132]68]. This dysregulated Wnt signaling can modulate the expression of the HOX family of proteins, including HOXB2 [[133]69]. GBM’s continuous modulation of intracellular and extracellular signals could lead to changes in HOXB2 expression, and this feedback mechanism helps explain the bidirectional causal relationship between HOXB2 and GBM. However, current research on HOXB2 and the ECM in the context of GBM is relatively limited, making it difficult to fully understand the deeper biological mechanisms involved. This represents a valuable research direction for future studies. Despite recent advancements in treatment approaches, the selection of therapeutic drugs for GBM remains relatively disappointing and challenging. The primary reason for this is the presence of the blood–brain barrier, which prevents most drugs from being accurately delivered to the target site. In our exploration of plasma proteins as potential pathogenic markers for GBM, although the evidence remains preliminary, our findings suggest that plasma proteins may be worth further investigation. The detection of proteins associated with GBM could present promising drug targets. Among the targets identified, GPX7 and CXCL10 emerged as potential candidates. Although there are no specific reports on the development of drugs targeting GPX7 for GBM, there is evidence suggesting that low-dose metformin upregulates ER-localized GPX7 to prevent cellular senescence [[134]70]. This highlights an interesting consideration: existing drugs that have already demonstrated efficacy may offer opportunities for application in different diseases, a question worth contemplating. Regarding CXCL10, while there are few reports on drug development targeting this protein, many studies suggest that CXCL10 could be a promising target in melanoma treatment. This is largely due to the dual role of the CXCL10/CXCR3 axis in promoting and counteracting cancer activity in various tissues and cells, particularly within the melanoma cells and their microenvironment [[135]71]. We believe that CXCL10 could be a promising drug target in glioma treatment, and our study provides preliminary evidence supporting this hypothesis. In future studies, we propose designing or screening small-molecule inhibitors or modulators targeting these two proteins. Virtual screening methods could be employed to identify known compounds capable of specifically binding to and inhibiting the activities of these proteins. Subsequent in vitro and in vivo experiments would help evaluate the effects of these compounds on GBM cell proliferation and migration, further assessing their potential clinical value. Given the role of CXCL10 in immune cell recruitment, it is likely associated with immune evasion mechanisms in GBM. This raises the possibility of developing combinational therapeutic strategies that integrate immune checkpoint inhibitors (e.g., PD-1/PD-L1 inhibitors) with CXCL10-related immune pathways to treat GBM. Recent studies have shown that a low-intensity focused ultrasound (LIFU)-guided sequential delivery strategy has been developed to enhance CD8 T cell infiltration and activity in GBM regions. This strategy involves the continuous delivery of CXCL10 to recruit CD8 T cells, along with interleukin 2 (IL-2) and aPD-L1 to reduce T cell exhaustion [[136]72].On the other hand, the antioxidant function of GPX7 might be linked to tumor cell resistance. Targeting GPX7 may help overcome GBM cell resistance to chemotherapy or radiotherapy. These biomarkers not only provide a means to evaluate drug efficacy but also offer a foundation for developing personalized therapeutic approaches. High-throughput screening technology to identify potential small molecules that interact with GPX7 from existing compound libraries is a potentially feasible approach. Once the compounds are screened, structure–activity relationship (SAR) analysis can be conducted to optimize the small molecules' structure, improving their selectivity and activity. Drug activity, solubility, and stability can be optimized through methods such as modifying functional groups or adjusting stereochemistry. Finally, in vitro and in vivo experiments can be performed to evaluate the potential therapeutic efficacy of the optimized compounds [[137]73]. It is worth noting that recent advancements in bioinformatics, including genomics and proteomics, have been complemented by the rapid development of technologies like machine learning and artificial intelligence. While MR analysis provides a straightforward method for causal inference, its capabilities are inherently limited. Incorporating other advanced techniques may improve the reliability of research findings. Recent studies highlight significant progress in applying deep learning techniques to medical imaging and genomics, demonstrating their potential for multidimensional data analysis. For example, an interpretable deep learning model emphasizing model transparency and interpretability has been validated in skin lesion classification [[138]74], providing valuable insights for analyzing GBM genomic and proteomic data. Applying such approaches to GBM research could enhance the precision of data interpretation and improve the transparency of causal relationship analyses. Additionally, a recent study introduced deep residual network technology [[139]75], which improves classification accuracy through multi-level feature extraction. This underscores the advantages of deep learning in processing high-dimensional data. Applying this approach to GBM gene and protein data analysis could help uncover deeper insights into disease mechanisms and optimize target screening efforts. Limitations MR studies typically require large sample sizes, and this issue is further compounded by the rarity of gliomas. To ensure that MR analyses provide robust references for drug target development, the current data