This study was granted an exemption from ethical review by the Medical Ethics Committee of The First People's Hospital of Zhaoqing (Reference No. B2026-08-03). The study used retrospectively collected, de-identified summary-level genetic and transcriptomic data, including Data Access Committee-approved controlled-access data from EGAD00010001657 and datasets obtained from GEO, eQTLGen, and GTEx in accordance with their applicable access and use conditions. No new participants were recruited, no biospecimens were collected, and no identifiable participant-level data were accessed. Ethical approval and informed consent for the original studies were obtained by the respective data generators, and the controlled-access data were used in accordance with the applicable Data Access Agreement.
Study design
This study employed a multistage framework to prioritize ferroptosis-related genes associated with glioblastoma (GBM) susceptibility and to evaluate their disease-relevant transcriptional effects (Figure 1). First, ferroptosis-related genes curated from FerrDb V2 were evaluated using two-sample Mendelian randomization (MR) with cis-expression quantitative trait locus (cis-eQTL) data and a large GBM genome-wide association study (GWAS). Inverse-variance weighted (IVW) MR served as the primary screening method, Bayesian weighted Mendelian randomization (BWMR) provided a complementary assessment of robustness, and an independent eQTL dataset was used for replication-stage MR. Second, genes supported by the genetic analyses were evaluated across three independent tumor–control transcriptomic cohorts, followed by cross-cohort meta-analysis. Third, patient-derived single-cell transcriptomic data were used to perform virtual gene perturbation in malignant cells and identify reproducible downstream transcriptional responses. These responses were subsequently characterized through functional enrichment and shared-network analyses. Overall, the genetic analyses were designed to identify genes associated with GBM susceptibility, whereas the transcriptomic and single-cell analyses assessed biological concordance and generated hypotheses for subsequent experimental validation.
Data sources
Ferroptosis-related genes were obtained from FerrDb V2, yielding 483 unique human protein-coding genes after gene-symbol harmonization and removal of duplicate entries14. Whole-blood cis-expression quantitative trait locus (cis-eQTL) summary statistics from the eQTLGen Consortium were used as the exposure dataset for discovery-stage Mendelian randomization (MR), whereas GTEx release V10 whole-blood cis-eQTL data served as an independent exposure dataset for replication-stage MR11,15. GBM outcome associations were obtained from controlled-access genome-wide association study (GWAS) summary statistics available through the European Genome-phenome Archive, comprising 6,183 cases and 18,169 controls of European ancestry16. Tissue-level expression of genetically prioritized genes was evaluated across three independent Gene Expression Omnibus (GEO) cohorts: GSE196533, comprising 61 grade 4 glioma samples annotated as GBM in the deposited metadata and nine non-neoplastic brain samples; GSE4290, comprising GBM and epilepsy-derived non-tumor brain samples; and GSE116520, containing paired tumor-core and peritumoral specimens together with non-neoplastic controls17,18,19. Patient-derived Smart-seq2 data from GSE131928 were used for malignant-cell virtual perturbation analysis20. Dataset characteristics and their respective analytical roles are summarized in Table 1. All analyses were performed using previously collected, de-identified datasets for which ethical approval and informed consent had been obtained in the original studies.
Genetic instrument selection and data harmonization
Candidate instruments were restricted to cis-expression quantitative trait loci (cis-eQTLs) associated with gene expression at genome-wide significance (P < 5 × 10⁻8). Variants were clumped using the European reference panel from the 1000 Genomes Project with a linkage disequilibrium (LD) threshold of r2 < 0.001 within a 10,000-kb window. Genes retaining fewer than three independent instruments after clumping were excluded from the primary Mendelian randomization (MR) analysis. Given the exploratory screening objective, a minimum of three instruments was prespecified to retain genes with sparse but strong cis-eQTL support while permitting multi-instrument IVW estimation. This threshold was accompanied by stringent LD clumping and F-statistic filtering; estimates based on only three or four instruments were interpreted cautiously, and sensitivity analyses were performed only when methodologically applicable. Instrument strength was assessed for each variant using the F statistic (F = β2/SE2), where β and SE represent the cis-eQTL effect estimate and its standard error, respectively. Variants with F < 10 were excluded to minimize weak-instrument bias21,22. Exposure and outcome datasets were harmonized by aligning effect alleles and effect directions. Duplicate variants, variants unavailable in the GBM GWAS dataset, and variants with incompatible allele coding were excluded. Because effect-allele frequencies were unavailable for the GBM GWAS, palindromic variants with ambiguous strand orientation were removed rather than inferred. For the same reason, a formal Steiger directionality test was not performed.
Mendelian randomization analyses
The association between genetically predicted gene expression and GBM susceptibility was evaluated using two-sample Mendelian randomization (MR). In the discovery stage, only genes with at least three independent cis-expression quantitative trait locus (cis-eQTL) instruments were included, and the inverse-variance weighted (IVW) method served as the primary analytical approach. Effect estimates were reported as odds ratios (ORs) with 95% confidence intervals (CIs) per unit increase in genetically predicted gene expression. To account for multiple testing across the genes evaluated, IVW P values were adjusted using the Benjamini–Hochberg procedure23. Genes with P < 0.05 and a Benjamini–Hochberg false discovery rate (BH-FDR) < 0.20 were retained as exploratory candidates. This relatively permissive FDR threshold was selected to reduce premature exclusion of potentially relevant genes during the discovery stage; therefore, candidate status was interpreted in conjunction with subsequent analyses rather than as confirmatory evidence. Bayesian-weighted Mendelian randomization (BWMR) was subsequently applied to the discovery-stage candidates using the same harmonized instruments24. Concordance between the IVW and BWMR results was assessed based on both statistical significance and the direction of the effect. Where permitted by the number of available instruments, Cochran's Q test, the MR-Egger intercept test, MR-PRESSO, and leave-one-out analyses were performed to evaluate heterogeneity, horizontal pleiotropy, influential outlying variants, and the influence of individual single-nucleotide polymorphisms (SNPs)22.
Replication-stage MR was performed using GTEx V10 whole-blood cis-eQTL data. The Wald ratio method was applied when only a single instrument was available, whereas the IVW method was used for genes with two or more instruments. Replication was defined by P < 0.05, BH-FDR < 0.20, and an effect direction consistent with the corresponding discovery-stage estimate. Because several GTEx replication estimates were based on only one or two instruments, they were interpreted as supportive evidence for replication rather than as independent evidence of causality.
Cross-cohort transcriptomic evaluation
The four genes prioritized by the discovery- and replication-stage Mendelian randomization (MR) analyses were evaluated across three independent transcriptomic cohorts. For GSE196533, raw RNA-sequencing count data were analyzed using DESeq225. Genes with counts below 10 in all but two samples were excluded, whereas the four target genes were retained regardless of expression filtering. Differential expression was assessed between 61 grade 4 glioma samples annotated as GBM in the deposited metadata and nine non-neoplastic brain samples.
For GSE4290, four samples lacking an explicit histopathological diagnosis were excluded, leaving 77 GBM and 23 non-tumor brain samples. Processed microarray intensities were log2-transformed, quantile-normalized, and analyzed using robust empirical Bayes linear models implemented in limma26. When multiple probes mapped to the same gene, the probe with the highest mean expression across all included samples was selected independently of differential expression significance.
GSE116520 comprised paired tumor-core and peritumoral specimens from 17 patients together with eight non-neoplastic controls. The deposited log-transformed, quantile-normalized expression data were analyzed using limma. Within-patient correlations between tumor-core and peritumoral samples were accounted for using patient-level blocking and the duplicateCorrelation function. Tumor-core versus control was the prespecified comparison for the cross-cohort meta-analysis, whereas peritumoral comparisons and the ordered control–peritumoral–tumor-core trend were evaluated separately.
Study-specific log2 fold changes and standard errors were pooled using a restricted maximum-likelihood random-effects model with Hartung–Knapp inference, as implemented in metafor. Between-study heterogeneity was assessed using Cochran's Q statistic and I2. Pooled P values for the four target genes were adjusted using the Benjamini–Hochberg procedure. Strong transcriptomic support was defined by a meta-analysis false discovery rate (FDR) < 0.05, cohort-level FDR significance in at least two datasets, and concordant effect directions across all three cohorts.
Single-cell virtual knockout analysis
Patient-derived Smart-seq2 data from GSE131928 were used to evaluate the four MR-replicated genes in malignant cells. Adult malignant cells were identified according to the annotations provided in the original study, and equal numbers of cells were randomly sampled from each eligible patient to minimize imbalance in patient representation. Virtual knockout was performed using scTenifoldKnk and repeated across five independent runs. Genes with a Benjamini–Hochberg-adjusted P < 0.05 in an individual run were considered significant. Downstream genes reproduced in at least three of the five runs were defined as the primary consensus set, whereas a more stringent four-of-five criterion was used for sensitivity analysis. These results were interpreted as computational predictions of transcriptional perturbation rather than evidence of direct molecular regulation.
Functional enrichment and shared-network analysis
Functional enrichment analysis was performed using target-specific downstream genes reproducibly identified in at least three of five virtual knockout runs. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathway enrichment were assessed using one-sided hypergeometric tests, with the 1,004 genes included in the network inference serving as the background gene set. P values were adjusted separately for each annotation database using the Benjamini–Hochberg procedure, and an adjusted P < 0.05 was considered statistically significant.
A bipartite network was constructed to represent the relationships between the four knockout targets and their consensus downstream genes. Genes associated with multiple targets were identified based on their shared degree, and overlap among target-specific gene sets was quantified using intersection counts and Jaccard indices. Network edges represent associations between reproducible computational perturbations and should not be interpreted as evidence of direct molecular interactions.
Statistical analysis and reproducibility
Unless otherwise specified, statistical tests were two-sided, and multiple comparisons were controlled using the Benjamini–Hochberg procedure. Analysis-specific significance criteria are described in the corresponding subsections. All analyses were performed using R or Python. Randomized procedures used prespecified seeds, and the analytical code, software versions, and detailed parameter settings were archived to support reproducibility. All datasets were previously collected and de-identified; ethical approval and informed consent were obtained in the original studies.