The dataset comprises a total of 176 GBM samples, categorized into 5 non-cancerous adjacent tissues and 171 tumor tissues, along with their respective clinical details from TCGA. Additionally, we obtained 59 GBM samples with matched clinical data from the Gene Expression Omnibus (GEO) repository (accession number GSE4421, https://www.ncbi.nlm.nih.gov/geo/), which were also subjected to normalization procedures. eQTL data encompassing 19,942 genes (refer to Supplementary Table 1), where significant single nucleotide polymorphisms (SNPs) were identified based on a significance threshold of p 5 $\times$ 10^-6. Furthermore, GWAS data for GBM from the FinnGen R10 release (https://www.finngen.fi/fi) were utilized, consisting of 253 GBM cases and 314,193 control subjects.

To guarantee the reliability of the MR analysis, potential instrumental variables (IVs) were chosen using a significance criterion of p 5 $\times$ 10^-6. This threshold was chosen due to the relatively small sample size of eQTL data (less than 50,000 samples per eQTL), which warranted a more stringent filtering criterion of p < 5 $\times$ 10^-6 to identify strong genetic instruments. To reduce the impact of linkage disequilibrium (LD) and ensure the independence of IVs, using parameters of R² = 0.001 and kb = 10,000. Additionally, palindromic SNPs, which could cause ambiguity in allele orientation, were excluded to avoid confounding results. Heterogeneous eQTLs were excluded. Pleiotropy in eQTLs was tested using the “mr_pleiotropy_test” function, and pleiotropic eQTLs were removed. Finally, SNP outliers were detected using the “run_mr_presso” function, and biased SNPs were excluded.

In the MR analysis, the core assumptions were carefully considered: the relevance of the IVs to the exposure, their independence from confounding factors, and the exclusion of direct effects on the outcome. The Wald ratio method was employed to obtain unconfounded MR estimates for each SNP. To further assess the robustness of the findings and account for potential pleiotropy, meta-analysis of the individual estimates was conducted using multiple methods. These methods allow for the assessment of causality while minimizing bias from pleiotropic effects and ensuring the reliability of the MR estimates.

Functional enrichment of differentially expressed gene modules between TCGA-derived high- and low-risk groups was performed using Gene Ontology (GO, https://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/kegg/) analyses via the clusterProfiler (v4.10.0, https://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) package in R. GO terms were identified with the enrichGO function (q 0.05), and KEGG pathways with the enrichKEGG function.

Expression data of intersecting genes were extracted from TCGA and combined with clinical data. The coxph function to select candidate genes significantly affecting GBM patient prognosis. LASSO regression analysis was then performed using the “glmnet” and “survminer” packages (v4.1-8, v0.4.9) in R language (v4.3.2) with the glmnet and cv.glmnet functions. Multivariate Cox regression analysis determined the final prognostic genes (p 0.05), and riskScore were calculated using the coef function (riskScore = $\mathrm{\Sigma }$ (coefficient_i $\times$ expression of signature gene_i)) (Supplementary Table 2). GBM samples were generated with the survival, regplot, and rms packages (v6.8-2, v1.1, v6.8-2) in combination with the coxph function. Differential expression of candidate genes was assessed.

All configured with the same parameter settings from a unified configuration file. The ssGSEA scores for each pathway were calculated using the gsva function and subsequently normalized via the normalize function.

GBM samples from TCGA with higher scores indicating stronger immune escape capability and poorer immunotherapy efficacy. Since the TIDE database (http://tide.dfci.harvard.edu/) does not include a specific category for GBM, the tumor type was designated as “Other” and the prior immunotherapy type was specified as “None”. Prior to uploading to the TIDE database, the TCGA GBM expression data were log2-transformed. The chi-square test was utilized to assess differences in immunotherapy response rates between these two groups.

The “pRRophetic” R package (v0.5) was employed to forecast sensitivity to chemotherapy drugs, and a correlation analysis was carried out between the drug sensitivity data and the riskScore, applying a filter with a significance threshold of p 0.001.

RSL3, ferrostatin-1, deferoxamine, JKE-1674, and atezolizumab (HY-100218C, HY-100579, HY-B1625R, HY-138153, HY-P9904) were purchased from MedChemExpress (Shanghai, China). The siRNAs and shRNA plasmids for DCTD and RRAS were constructed by Shanghai Genechem (China). Antibodies against RRAS (27457-1-AP), SLC7A11 (26864-1-AP), DCTD (68357-1-Ig), ACSL4 (22401-1-AP), GAPDH (10494-1-AP), and GPX4 (30388-1-AP) were obtained from Proteintech (Wuhan, China). All the primary antibodies were diluted at a ratio of 1:1000. Specifically, the sequence for both si-RRAS1 and sh-RRAS1 is 5^′-ACACGAAGATCTGCAGTGTGGAT-3^′, while the sequence for si-RRAS2 and sh-RRAS2 is 5^′-AGACGAAGATCTGCAGTGTGGAC-3^′. Additionally, the sequences for si-DCTD1 and sh-DCTD1 are 5^′-GGAUCUUCUGCUUCCAAUA-3^′, and those for si-DCTD2 and sh-DCTD2 are 5^′-CCAAGACAGAUCCUGUUAU-3^′.

Human glioblastoma cell lines U-87MG, A172, and TG-905 (ATCC, Shanghai, China) were cultured in a CO₂ incubator (TCI-80, Panasonic Healthcare Holdings Co., Ltd, Ora, Japan) at 37 °C and monitored every 8 or 24 h as required. Cells were subcultured at 80–95% confluence. All disposable materials were UV-sterilized in a laminar flow cabinet (Jiesen, Xuzhou, Jiangsu, China), and procedures were performed with sterile gloves. After trypsin digestion, cells were suspended in complete growth medium, centrifuged, and resuspended in fresh medium for continued culture. During cell culture, we regularly add mycoplasma scavengers to the medium to avoid mycoplasma contamination.

Cells were cultured in 6-well plates to reach 70–80% confluence. The medium was replaced with serum-free opti-MEM before transfection. Transfection was performed using Lipofectamine 3000 (L3000015, Thermo Fisher Scientific, Shanghai, China) according to the manufacturer’s instructions, with either plasmid DNA or siRNA. This complex was then added dropwise to the cells in the 6-well plates. For siRNA transfection, a similar procedure was followed, but with siRNA molecules instead of plasmids. The cells were exposed to the transfection complexes for a 6-hour period at 37 °C. Following this incubation, the existing medium was substituted with new complete medium supplemented with 10% fetal bovine serum (FBS), aiming to promote cell recovery. Post-transfection, cells were incubated for 48 hours under standard culture conditions before being used for further experimentation, such as RNA/protein extraction, phenotypic analysis, or functional assays. All cell lines were validated by short tandem repeat (STR) profiling and tested negative for mycoplasma.

Protein extracts were prepared by lysing cells in cold lysis buffer with protease and phosphatase inhibitors, followed by centrifugation to obtain the soluble protein fraction. Protein concentrations were determined using a BCA assay (E-BC-K318-M, Elabscience Biotechnology Inc., Wuhan, Hubei, China), and samples were mixed with sodium dodecyl sulfate (SDS) buffer and denatured by heating. Proteins were separated by sodium dodecyl sulfate - polyacrylamide gel electrophoresis (SDS-PAGE) using a 10% resolving gel.

Cell viability was evaluated using the CellTiter-Glo 3D assay by measuring luminescence after treatment with ferroptosis inducers and inhibitors. Lipid peroxidation was assessed using the BODIPY C11 probe, with the ratio of oxidized to reduced fluorescence indicating peroxidation levels. Glioblastoma cells were treated for 24–72 hours. Luminescence and fluorescence were measured using microplate readers. Higher luminescence indicates higher viability, while a higher oxidized/reduced fluorescence ratio indicates greater lipid peroxidation.

Regarding the animal housing environment, experimental animals are housed in a dedicated animal facility with stringent environmental controls. This facility complies with the CV grade standard. The indoor temperature is consistently maintained at 22 $\pm$ 2 °C, a range that has been extensively validated. This temperature range effectively minimizes physiological variations caused by fluctuations, thereby ensuring the stability and reliability of experimental results. Humidity is controlled at 50% $\pm$ 5%, which helps prevent respiratory diseases and ensures animal health. The lighting environment was kept on a 12-hour light/dark cycle to simulate natural circadian rhythms and maintain the animals’ biological regularity. In the xenograft experiments, 1 $\times$ 10⁶ U-87MG glioblastoma cells, transfected with siRNA or control siRNA (sicon), were mixed with 100 µL Matrigel (1:1) and implanted subcutaneously into the right flank of 5–6-week-old male nude mice. Tumor growth was monitored every 2–3 days by caliper measurements and volume calculated using V = L $\times$ W² $\times$ 0.52. Mice were assigned to treatment groups once tumors reached 80–100 mm³. They were administered JKE-1674 (25 mg/kg) orally every other day, and atezolizumab (10 mg/kg) intravenously twice weekly, for a duration of 3 weeks. Throughout the treatment period, peripheral blood was collected weekly for monitoring liver function, specifically alanine transaminase (ALT) and aspartate aminotransferase (AST) levels, to assess potential drug toxicity. Following the completion of the 3-week treatment phase, the mice were euthanized via cervical dislocation, and their tumor tissues were collected for subsequent analysis, including histological evaluation, protein analysis, and RNA extraction. Our in vivo studies strictly follow the “ARRIVE and PREPARE” guidelines, which promise to improve the scientific nature and transparency of the study while ensuring ethical treatment of the animals.

For normally distributed variables, an independent Student’s t-test was used to assess statistical significance, whereas the Mann-Whitney U test (also known as the Wilcoxon rank-sum test) was used for non-normally distributed data.

We analyzed eQTL data from blood samples for 19,942 genes to explore associations with GBM. After filtering, 15,306 eQTLs and 52,674 SNPs were selected for MR analysis. Cross-referencing with GBM GWAS data identified 49,823 SNPs with F-statistics $>$10. A total of 15,035 genes were included in MR analysis. Using inverse-variance weighted (IVW) method, 250 genes showed significant associations with GBM risk (147 increased risk, 103 reduced risk) (Table 1, Supplementary Tables 3–7).

To further screen for GBM-related genes, a differential analysis was conducted on mRNA data derived from TCGA cancerous tissues and adjacent normal tissues. Using a false discovery rate (FDR) threshold of less than 0.05 and a logFC value greater than 0.585 as screening criteria, a total of 9942 differentially expressed genes were identified (Supplementary Table 8). Heatmap and volcano plot analyses revealed significant differences in gene expression between the two groups, with fewer upregulated genes observed in the tumor group compared to the adjacent normal group (Fig. 1A,B). By intersecting the risk gene sets obtained from both TCGA and MR analyses, we identified 50 genes that were significantly associated with GBM. Among these, 39 genes were classified as high-risk genes for GBM, while 11 genes were identified as protective genes for GBM (Fig. 1C,D, Supplementary Table 9).

Fig. 1.

Difference analysis of GBM mRNA expression between cancer group and paracancer group, and acquisition of risk genes. (A,B) Heatmap and volcano plot illustrating genes differentially expressed between the cancer and paracancer groups. Elevated expression in tumor tissues is indicated in red, while higher expression in normal tissues is shown in blue. (C) Overlap between genes highly expressed in the cancer group and those with an odds ratio (OR) greater than 1 based on MR analysis. (D) Overlap between genes predominantly expressed in the paracancer group and those with an odds ratio (OR) less than 1 according to MR analysis. MR, Mendelian randomization; GBM, glioblastoma; DEG, differentially expressed genes.

Univariate Cox regression analysis conducted on the 50 candidate genes indicated that 11 genes exhibited a statistically significant link to the prognosis of GBM patients (Fig. 2A). Subsequent MR analysis indicated that the eQTLs of these 11 genes may have a significant impact on GBM prognosis (Fig. 2B). The Circos plot visually illustrates the chromosomal locations of these 11 genes (Fig. 2C). Cross-validation curve analysis and LASSO coefficient path analysis determined four risk genes (FSTL1, FXYD5, RRAS, and RNF216P1) as variables in the model and provided the model coefficients for each gene (Fig. 2D,E, Supplementary Table 10).

Fig. 2.

Acquisition of GBM candidate genes and establishment of LASSO regression model. (A) Univariate Cox analysis of candidate genes in GBM patients (p 0.05). (B) eQTL MR Analysis in patients with GBM, this figure shows the OR values of SNP and GBM for eQTL genes in IVW and Weighted median methods. (C) The chromosomal circle shows the location of 11 risk genes. (D,E) LASSO regression analysis cross-validation curves and LASSO coefficient path diagram in LASSO regression analysis. eQTL, expression quantitative trait loci; IVW; inverse-variance weighted.

TCGA tumor samples were divided randomly into training and testing cohorts. Four key genes were significantly higher in the high-risk group than in the low-risk group. High-risk individuals showed higher mortality and shorter survival times (Fig. 3A–C). Survival analysis confirmed that high-risk patients had notably lower survival rates (Fig. 3D–F). Data validation analysis was also performed using the GSE4421 dataset. Unfortunately, the RNA expression data in this dataset does not include RNF216P1, so the results from this analysis cannot be used (Supplementary Fig. 1). These findings were consistent when validated across all TCGA tumor samples, underscoring the importance of these four genes in GBM prognosis.

Fig. 3.

Survival status and gene expression differences across high- and low-risk categories. ((A) Training cohort, (B) Testing cohort, (C) Combined cohort) The expression levels of four critical genes and their association with riskScore and survival outcomes. ((D) Training cohort, (E) Testing cohort, (F) Combined cohort) Survival regression analysis comparing high- and low-risk groups.

In conjunction with age and gender, was significantly associated with GBM prognosis, exhibiting a hazard ratio (HR) greater than 1.5. This suggests that an elevated riskScore is correlated with a poorer prognosis (Fig. 4A,B). Receiver operating characteristic (ROC) curves and the concordance index (C-index) further demonstrated that the riskScore possessed high specificity and sensitivity for predicting GBM prognosis (Fig. 4C,D, Supplementary Fig. 2). A more sophisticated predictive model integrating age, gender, and the riskScore was subsequently developed to enhance prediction accuracy (Fig. 4E). This model exhibited robust accuracy in forecasting 2- and 3-year survival rates, albeit with slightly diminished performance for 1-year survival prediction (Fig. 4F). Principal component analysis (PCA) effectively differentiated between high-risk and low-risk groups, particularly through the application of the four identified risk-associated genes (Fig. 4G–I).

Fig. 4.

Validation of risk regression model. (A) Univariate regression analysis was conducted to assess the impact of riskScore, sex, and age on glioblastoma (GBM) patients. (B) A multivariate regression model was employed to evaluate the combined effects of riskScore, gender, and age in GBM patients. (C) Receiver operating characteristic (ROC) curves were generated for riskScore, sex, and age. (D) ROC curves specifically for riskScore at 1, 2, and 3 years were plotted. (E) A nomogram incorporating gender, age, and riskScore was developed for GBM patients. *p 0.05; ***p 0.001. (F) Calibration curves were constructed to assess the predictive accuracy of the model for GBM patients at 1, 2, and 3 years. (G) Principal component analysis (PCA) was performed to compare the gene expression profiles of the high-risk and low-risk groups. (H) PCA was also applied to analyze the fifty overlapping risk genes between the two groups. (I) Additionally, PCA was utilized to evaluate the expression patterns of four key gene clusters in the high-risk and low-risk cohorts.

These genes were categorized into two groups: those more highly expressed in the high-risk group and those more prevalent in the low-risk group. GO and KEGG pathway enrichment analyses were then conducted for each gene group. Results revealed in the high-risk group, that gene sets, including those linked to cytokine activity, chemokine activity, and immune cell migration were significantly enriched. Conversely, in the low-risk group, gene sets associated with ion transport and neurotransmitter regulation, such as modulation of chemical synaptic transmission and regulation of metal ion transport, were mainly enriched (Fig. 5A–D, Supplementary Tables 11–13). Consistent with these findings, KEGG enrichment analysis mirrored the GO results, demonstrating high-risk group was related to cytokines, chemokines, and inflammation, such as interleukin-17 (IL-17) and tumor necrosis factor (TNF) signaling (Fig. 5E,F, Supplementary Table 14). The low-risk group was significantly associated with neuroactive ligand-receptor interactions and calcium signaling pathways, which are involved in ion transport and neurotransmitter regulation. The low-risk group exhibited increased expression of genes related to ion transport and neurotransmitter regulation (Fig. 5G,H, Supplementary Table 15). These findings indicate that the riskScore effectively distinguishes major gene expression patterns in GBM patients.

Fig. 5.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses in high-low risk groups. (A,B) GO analysis in high-risk groups. (C,D) GO analysis in low-risk groups. (E,F) KEGG analysis in high-risk groups. (G,H) KEGG analysis in low-risk groups. Sector diagram: Different sector blocks represent different enrichment pathways, and the height of the sector in the inner circle represents the proportion of gene sets where the enriched genes are located. The higher the sector is, the higher the proportion is.

Next, we conducted a comprehensive examination of genetic mutation characteristics in both high- and low-risk groups based on TCGA data. The findings revealed notable variations in the mutation rates of the wild-type human p53 (TP53), epidermal growth factor receptor (EGFR), and phosphatase and tensin homolog (PTEN) across the two groups. In particular, in the high-risk group, mutations in PTEN were more prevalent, but mutation frequencies of TP53 and EGFR were more prevalent in low-risk group (Fig. 6A,B). Co-mutation analysis further revealed notable patterns of co-occurrence and mutual exclusivity. For instance, PTEN and TP53 co-occurred more frequently in the high-risk group, while TP53 co-occurred with isocitrate dehydrogenase 1 (IDH1) and alpha-thalassemia mental retardation X-linked (ATRX) in both groups. Conversely, TP53 and EGFR exhibited mutual exclusivity. Notably, the high-risk group exhibited a larger number of concurrently occurring gene pairs (16 pairs) in comparison to the low-risk group (9 pairs) (Fig. 6C,D). variant allele frequency (VAF) analysis revealed that the high-risk group had generally lower average VAF values than the low-risk group. Among the high-risk group, PTEN, TP53, and erythrocytic 1 (SPTA1) were the top genes based on VAF. In contrast, the low-risk group showed ATRX, retinoblastoma protein 1 (RB1), and TP53 as the most prominent genes. Importantly, the average VAF of TP53 differed significantly between the two groups (Fig. 6E,F).

Fig. 6.

Gene mutation characteristics and tumor mutational burden (TMB) status of high and low-risk groups. ((A) High-risk cohort, (B) Low-risk cohort) The waterfall plot illustrates the mutation frequency and types of genes within the two groups. The upper section of the plot represents TMB, while the lower section provides clinical details and risk classifications of the patients. ((C) High-risk cohort, (D) Low-risk cohort) The co-mutation patterns and mutation frequencies of the top 20 genes in both groups are displayed. Green highlights co-occurring mutations, whereas yellow signifies mutually exclusive mutations. ((E) High-risk cohort, (F) Low-risk cohort) Box plots depicting the mean VAFs observed in the two groups. (G) A violin plot demonstrating the correlation between the high-risk cohort and TMB levels. (H) Survival analysis comparing groups with high and low TMB. (I) Survival curves stratifying patients based on both TMB levels and riskScore categories. VAF, variant allele frequency.

We further explored the impact of tumor mutational burden (TMB) on GBM and found that patients with higher TMB levels had notably better survival outcomes than those with lower TMB levels (Fig. 6H). Nevertheless, no significant difference in TMB was observed (Fig. 6G). The patients were divided into four groups based on TMB levels and riskScores. The low TMB but high-risk group had the worst survival prognosis. The high TMB and low-risk group showed the better quality of life (Fig. 6I). Collectively, these findings indicate that the riskScore effectively distinguishes gene mutation patterns and characteristics among the different groups.

We utilized the R programming language to assess immune infiltration levels in the samples through multiple immune cell estimation approaches, such as Estimate, ssGSEA, MCPcounter, Timer, Epic, CIBERSORT, and Quantiseq. Out of these seven techniques, the high-risk group demonstrated markedly increased infiltration levels of T cells (specifically CD4+ and CD8+ T cells) when compared to the low-risk group. Furthermore, in nearly all methods except CIBERSORT, the high-risk group displayed significantly greater infiltration of monocytes and tumor-associated fibroblasts. (Fig. 7A–G, Supplementary Table 16). The high-risk group also exhibited higher stromal and immune scores overall. Additionally, scores associated with immune inflammatory responses—such as those involving major histocompatibility complex (MHC), interferon (IFN), antigen-presenting cell (APC), and T cell-related gene sets—were considerably elevated in this group (Fig. 7D,H, Supplementary Table 16). Using the ssGSEA method, the infiltration levels of myeloid-derived suppressor cells (MDSCs) were found to be significantly higher in the high-risk group compared to the low-risk group (Fig. 7G, Supplementary Table 16). Examination of immune checkpoint expression revealed that several key checkpoints, such as cytotoxic T lymphocyte-associated antigen-4 (CTLA4), Programmed cell death protein 1‌ (PDCD1), Programmed Cell Death 1 Ligand 1 (CD274), and Recombinant Hepatitis A Virus Cellular Receptor 2 (HAVCR2), were significantly overexpressed in the high-risk group (Fig. 7I,J). Collectively, these results highlight a clear distinction in the tumor immune microenvironment between the high-risk and low-risk groups.

Fig. 7.

Relationship between high-low risk groups and tumor immune microenvironment. (A–G) The relationship between high and low-risk groups and the level of immune cell infiltration in 7 immune scoring tools. The vertical coordinates indicate different methods, and the horizontal coordinates indicate different immune cells. (H) Heatmap show the score levels of immune-related gene sets in all samples from both groups. The radar map showed the difference in the expression levels of co-inhibitory (I) and co-stimulatory immune checkpoints (J) in high and low-risk groups. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; *****p 0.00001; ns $>$ 0.05.

Study used the TIDE platform to judge immunotherapy response immune cell infiltration in GBM samples. Higher infiltration of CAFs and increased expression of IFN, CD8, and programmed cell death ligand 1 (PD-L1) in high-risk group. The high-risk group also had higher immune dysfunction scores, indicating more significant infiltration of tumor-associated immune cells. Unlike previous findings, the high-risk group had lower infiltration of MDSCs and M2 macrophages. However, there was no significant difference in predicted immunotherapy response between the two groups (Supplementary Fig. 3, Supplementary Table 17).

The study used the R package “pRRophetic” to estimate TCGA GBM patients’ sensitivity to various chemotherapeutic agents. A significant inversion was found in dasatinib, lapatinib, and temsirolimus, suggesting higher riskScores correlate with better responses to these drugs. Conversely, a positive correlation was observed between riskScore and the IC₅₀ values of elacridar, lestaurtinib, and tubastatin A, indicating lower riskScores may lead to greater sensitivity to these drugs (Supplementary Fig. 4). A systematic analysis of the link between riskScore and chemotherapy drug sensitivity allowed us to identify differing drug sensitivity among patient groups with varying riskScores.

The boxplots of gene expression demonstrated that the mRNA expression levels of the four risk-associated genes were significantly elevated in tumor tissues relative to adjacent normal tissues (Fig. 8A). Subsequently, we investigated the correlation between the abundance of tumor-infiltrating immune cells and the expression of these four risk genes to elucidate their roles in immune cell infiltration. Across six distinct immune cell infiltration scoring methodologies, FXYD5 exhibited a significant positive correlation with macrophage infiltration levels. Furthermore, FXYD5 demonstrated a positive relationship CD4+ and CD8+ T cells across the majority of the scoring methods employed. In the Epic and MCPcounter models, FXYD5 was found to be positively linked with CAFs. Additionally, using the ssGSEA, MCPcounter, and CIBERSORT approaches, a positive association was observed between FXYD5 and B cell infiltration levels (Fig. 8B–H). In a similar manner, RRAS expression levels were related to the presence of CAFs, T cells, and macrophages, within the microenvironment. However, RNF216P1 expression showed a bad relationship with multiple tumor-infiltrating immune cells, such as macrophages, CAFs, and T cells, as indicated by six different scoring systems (Fig. 8B–H). These results suggest that riskScore serves as a significant indicator for evaluating immune cell infiltration in GBM patients.

Fig. 8.

Relationship between key genes and tumor immune microenvironment. (A) Expression levels of 4 key genes in cancer and para-cancer tissues in TCGA GBM. ((B) Quantiseq, (C) Epic, (D) Timer, (E) Estimate, (F) MCPcounter, (G) CIBERSORT, (H) ssGSEA) Relationship between immune cell infiltration level and key genes in 7 kinds of immune score calculation tools. The orange line indicates a positive correlation between the gene and the level of immune cell infiltration, whereas the green line indicates a negative correlation. ***p 0.001.

The GBM samples were categorized into two groups—low expression and high expression—according to the median expression values of the four risk-associated genes, after which survival analysis was conducted. The findings revealed that higher expression levels of DCTD and RRAS were linked to worse prognostic outcomes (Fig. 9A,B), whereas the remaining two genes did not exhibit a statistically significant difference between the two expression groups (Fig. 9C,D). Based on these results, we selected DCTD and RRAS for further experiments. We transfected GBM U-87MG cells with siRNAs targeting DCTD and RRAS.

Fig. 9.

Survival analysis of four risk genes. (A–D) The KM survival curves of four survival genes. Based on these results, we selected DCTD and RRAS for further experiments. RRAS, Ras-related protein.

We transfected GBM U-87MG cells with siRNAs targeting DCTD and RRAS. Knockdown of DCTD or RRAS led to a reduction in their respective protein levels (Fig. 10A,B). MTT, colony formation, wound healing, and invasion assays demonstrated that knockdown of DCTD or RRAS reduced U-87MG cell proliferation, migration, and invasion (Fig. 10C–J, Supplementary Fig. 5). Using a CDX model, we found that tumors with DCTD or RRAS knockdown were significantly smaller than those in the control group (Fig. 10K–N), We have also added a new cell line, U251, and obtained similar results (Supplementary Fig. 6) indicating that knockdown of DCTD or RRAS inhibits GBM CDX model tumor growth.

Fig. 10.

DCTD and RRAS promote the proliferation, migration and invasion of GBM. (A,B) Protein expression after DCTD or RRAS knockdown. (C,D) MTT experiment after DCTD or RRAS knockdown. (E,F) Clone formation experiment after DCTD or RRAS knockdown. (G,H) Wound healing test after DCTD or RRAS knockdown. Scale bar = 200 µm. (I,J) Invasion Assay after DCTD or RRAS knockdown. Scale bar = 100 µm. (K–N) CDX model construction and tumor weight detection after DCTD or RRAS knockdown. (Related experiments were independently replicated three times.) **p 0.01; ***p 0.001; ****p 0.0001.

Using GEPIA for correlation analysis, a significant positive correlation was found between RRAS and GPX4, a key ferroptosis gene (Supplementary Fig. 7A). Although the correlation between DCTD and GPX4 was not statistically significant, the p-value was close to 0.05 with a positive correlation coefficient (Supplementary Fig. 7B). Knockdown of DCTD or RRAS led to decreased expression of ferroptosis inhibitory proteins SLC7A11 and GPX4, and increased expression of the ferroptosis-promoting protein ACSL4 (Fig. 11A–D, Supplementary Fig. 5A–D). Colony formation, wound healing, and invasion assays showed that cells with si-DCTD or si-RRAS and exogenous ACSL4 had lower proliferation, migration, and invasion capabilities compared to cells with si-DCTD or si-RRAS alone. Conversely, overexpression of DCTD or RRAS enhanced these abilities, and knocking down ACSL4 in DCTD or RRAS overexpression cells increased proliferation, migration, and invasion (Fig. 11E–P, Supplementary Fig. 5G–I). This indicates that DCTD and RRAS influence the ferroptosis signaling pathway and affect U-87MG cell phenotype. Meanwhile, we investigated the mechanism by which DCTD regulates ferroptosis. Initially, using q-PCR, we demonstrated that DCTD knockout did not induce changes in the transcriptional level of ACSL4 (Supplementary Fig. 8A). Subsequently, we hypothesized whether DCTD modulates ACSL4 expression by influencing ubiquitination-related enzymes. To this end, we predicted E3 ubiquitin ligases regulating ACSL4 via the ubibrowser website and identified autocrine motility factor receptor (AMFR) as a regulator of ACSL4 (Supplementary Fig. 8B). Furthermore, through analysis on gepia2.cancer, we observed a positive correlation between DCTD and AMFR (Supplementary Fig. 8C). Next, we transfected si-RRAS and si-DCTD into the U251 cell line and assessed the protein expression levels of RRAS, DCTD, and AMFR using Western blotting. Our results revealed that RRAS and DCTD regulate the protein level of AMFR, with a positive correlation among them (Supplementary Fig. 8D,E). Finally, we designed an experiment to transfect an AMFR overexpression plasmid following DCTD knockout and examined the ubiquitination level of ACSL4. We found that the ubiquitination level of ACSL4 decreased significantly upon DCTD knockout. However, when the AMFR overexpression plasmid was transfected subsequent to DCTD knockout, the ubiquitination level of ACSL4 was restored (Supplementary Fig. 8F).

Fig. 11.

Detection of ferroptosis and effect on GBM cell phenotype after knockdown or overexpression of DCTD and RRAS. (A,B) Expression of ferroptosis signaling protein after DCTD or RRAS knockdown. (C,D) Expression of ferroptosis signaling protein after DCTD or RRAS overexpression. (E–P) Clone formation experiment, migration experiment and invasion experiment of different groups of U-87MG. (I, K) Scale bar = 200 µm. (M, O) Scale bar = 100 µm. (Related experiments were independently replicated three times.) *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001.

To validate the effect of ferroptosis on GBM cells, we treated cells with the ferroptosis activator RSL3 and the ferroptosis inhibitors ferrostatin-1 and deferoxamine (iron chelator). Activation of the ferroptosis pathway significantly reduced cell viability in three GBM cell lines, and this effect was reversed by ferroptosis inhibitors (Fig. 12A–C), indicating that activation of the ferroptosis pathway decreases GBM cell viability. In U-87MG cells with DCTD or RRAS knockdown or overexpression, RSL3 treatment at varying concentrations for 24 hours showed that the IC₅₀ of RSL3 decreased with DCTD or RRAS knockdown and increased with overexpression. Cells treated with RSL3 showed reduced viability in si-DCTD, si-RRAS, OE-DCTD, and OE-RRAS groups compared to controls. Additionally, lipid peroxidation levels were higher in the DCTD or RRAS knockdown groups and lower in the overexpression groups after RSL3 treatment (Fig. 12D–O), suggesting that knockdown of DCTD or RRAS enhances ferroptosis in U-87MG cells treated with RSL3. We detected the intracellular reactive oxygen species (ROS) levels after the knockout and overexpression of RRAS and DCTD using a flow cytometer. We found that the expression level of ROS increased significantly after the knockout of RRAS and DCTD, while it decreased significantly after their overexpression (Supplementary Fig. 9).

Fig. 12.

Determination of the effect of DCTD and RRAS on the level of ferroptosis of GBM cells. ((A) TG-905, (B) U-87MG, (C) A172) Cell viability analyses of different groups of three GBM cells treated with ferroptosis activators or ferroptosis inhibitors. (D,E,G,H,J,K,M,N) Detection of IC₅₀ values of RAS-selective lethal small molecule 3 (RSL3) after knockdown or overexpression of DCTD, RRAS by cell viability analyses. (F,I,L,O) Detection of lipid peroxidation levels after application of RSL3 or Ferrostain-1 in knockdown or overexpression of DCTD and RRAS U-87MG cells. (Related experiments were independently replicated three times.) *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001.

Given that JKE-1647 is more easily metabolized in vivo and induces ferroptosis [18], we chose JKE-1647 as the ferroptosis inducer for in vivo experiments. Since DCTD and RRAS are positively correlated with immune cell infiltration, we selected the PD-L1 inhibitor atezolizumab to verify the efficacy of immunotherapy for GBM. We established a subcutaneous tumor model with three groups: si-con + JKE-1647, si-con + JKE-1647 + atezolizumab, and si-DCTD/si-RRAS + JKE-1647 + atezolizumab. Consistent with in vitro results, ferroptosis activation significantly inhibited subcutaneous tumor growth, and the combination of JKE-1647 and atezolizumab had a stronger inhibitory effect than JKE-1647 alone. As expected, knockdown of DCTD or RRAS combined with JKE-1647 and atezolizumab more effectively inhibited GBM subcutaneous tumor growth (Fig. 13A–C, F–H), indicating that DCTD or RRAS knockdown promotes ferroptosis activation. Plasma ALT and AST levels showed no significant differences among the groups, suggesting that JKE-1647 combined with atezolizumab does not cause noticeable hepatotoxicity (Fig. 13D,E,I,J). We added the murine glioma cell line GL261 to establish a subcutaneous tumor model, and obtained the same conclusion (Supplementary Fig. 10).

Fig. 13.

Effect of JKE-1647 combined with atezolizumab on the growth of subcutaneous transplanted tumors in nude mice. (A–C, F–H) CDX model construction and tumor volume and weight detection in different groups. (D,E,I,J) Detection of liver function in different groups of nude mice. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns $>$ 0.05.