Gene expression data and the glioma patient survival information in The Cancer Genome Atlas (TCGA) program were downloaded from the GlioVis GBM platform (http://gliovis.bioinfo.cnio.es/) (21), in which the RNA-seq data were processed by the normalized read counts from the preprocessed data with log2 transformation.

We downloaded 253 ferroptosis regulators, 111 ferroptosis markers, and 95 ferroptosis–disease associations from FerrDb (http://www.zhounan.org/ferrdb). FerrDb can be used to obtain insights into ferroptosis ( Figure 1A ) (22). The ferroptosis gene expression data in the TCGA-GBM were obtained, and the ferroptosis genes significantly correlated with prognosis were determined based on the Cox regression method (p-value < 0.05). The expression matrix of relevant iron metabolism genes was extracted. Then, the score of each TCGA-GBM sample was calculated according to the principal component analysis (PCA) algorithm. PCA was applied to establish a ferroptosis-relevant score model, namely, the PCA score, and the PCA score was determined based on the formula: PCA score = (PC1i + PC2i), in which i represents the level of ferroptosis family pattern-related signature genes. Both PC1 and PC2 were regarded as signature scores. According to the threshold, the patients were classified into high- and low-score groups, while the relevant survival curves were generated and checked by the log-rank test between the two groups.

We identified the differentially expressed genes (DEGs) between the high-score and low-score groups. The “limma” package (23) of R was used to identify DEGs. Fold-change (FC) >1 and adjusted p-value <0.05 were set as the judgment standard for DEGs. An adjusted p-value <0.05 was applied to control the false-positive rate (24).

Based on the screened DEGs, unsupervised clustering was applied to cluster distant genotypes of the TCGA-GBM samples by the “ConsensusClusterPlus” package (25). Heatmaps were drawn to visualize the clinical characteristics of different gene clusters. The GBM score of each individual was also determined by these ferroptosis-associated DEGs, which we called the PCA gene score. In terms of cutoff value, the cases were classified into the high-gene score and low-gene score groups. Kaplan–Meier (KM) curves were generated and compared with the log-rank test between the high-gene score and low-gene score groups.

Gene ontology (GO) analysis (26) was used to analyze the function of the gene products and to identify their biological feature. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (27) enrichment analysis was mainly used to explain gene function and was widely used in the bioinformatics research field. Based on DEG, it was mainly applied to the “ClusterProfiler” software (28) during the GO and KEGG analyses. The “gsva” package (29) was mainly used for single-sample gene set enrichment analysis (ssGSEA) to determine the proportion of infiltrating immune cells and also to determine the activity difference of immune pathways in different groups. The correlations of infiltrating immune cells and PCA score/gene score were shown in correlation heatmaps.

The milled tissue and cells were lysed and boiled in sodium dodecyl sulfate (SDS) sample loading buffer, separated by 10% SDS-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to PVDF membranes. The membranes were blocked in 5% skimmed milk for 1 h, washed three times with TBST, incubated with the appropriate primary antibody in TBST for 1 h, then washed twice, and incubated with horseradish peroxide-conjugated anti-immunoglobulin (1:5,000 dilution) for 1 h at room temperature. After three washes with TBST, the membranes were developed using ECL-enhanced chemiluminescence reagents.

Tissues were fixed in 4% paraformaldehyde for 24 h, followed by embedding in paraffin. Sections were cut to a thickness of 5 μm and incubated with anti-MAP1LC3A, anti-OLFM1, anti-CEND1, anti-CLTB, anti-PRKAR1B, anti-HRAS, anti-DDRGK1, anti-PITHD1, and anti-DEAF1 primary antibodies. Hematoxylin restaining was followed by separate photographic observation of each slide in a positive area of ×40 fields, which were independently analyzed by three previously uninformed pathologists.

To compare the overall survival (OS) of the subgroups, the KM curve with the log-rank test was applied (30). In order to determine the prognostic value of the iron production death gene in the research process, this paper conducted Cox regression analysis on the collected data and applied the Mann–Whitney method when comparing the degree of immune cell infiltration of the two groups. The correlations of infiltrating immune cells and PCA score/gene score were established using the Pearson correlation coefficient. The statistical analyses were carried out based on R (version 4.0.5).

Using univariate regression analysis, we screened 14 ferroptosis genes associated with OS in patients with GBM. According to the optimal threshold of the level of these genes in the TCGA-GBM, the cases were divided into two groups with high and low expression. The KM curve showed significant prognostic differences between the high and low gene expression groups (p-value < 0.05). The PCA algorithm was applied to determine the PCA score of each TCGA-GBM individual based on these 14 genes. The results showed that a high score has a positive relationship with prognosis (p = 0.002, Figure 1B ).

Up to 231 DEGs were identified, among which 183 were upregulated and 48 were downregulated in the high-score group. To further explore the biological role of these ferroptosis-related genes in GBM, we performed GO and KEGG analyses on these DEGs. The GO result reflected that the GO annotations of DEGs can be classified into three categories: biological processes (BPs), cell components (CCs), and molecular functions (MFs). The terms were arranged based on the false discovery rate (FDR) value, and further screening was performed when FDR <0.05. Screening showed that the DEGs were enriched in processes of regulation of long-term neuronal synaptic plasticity, proteasome-mediated ubiquitin-dependent protein degradation, vesicle-mediated transport in synapses, and proteasomal protein catabolic processes ( Figures S1A, B ). The signal pathway results were ranked based on FDR value and were screened under the condition of FDR <0.05. The KEGG result reflected that the DEGs were mainly enriched in endocytosis, Alzheimer’s disease, and pathways of neurodegeneration in multiple diseases ( Figures S1C, D ).

To study the association of the expression level and GBM subtypes of the 231 ferroptosis-related DEGs, we conducted univariate Cox regression and consensus clustering test with all GBM cases in the TCGA. When the clustering variable (k) increased from 2 to 10, it can be found that when k = 2, there appeared the highest intragroup associations and the lowest intergroup associations, which reflected that the 158 GBM cases can be perfectly divided into two groups by the 25 prognostic DEGs ( Figures 2A–D ). The clinical features, such as PCA score (high score or low score), sex (male or female), age (≤65 or >65 years), and survival condition, were listed in a heatmap ( Figure 3A ), and a few differences in sex and age between the two groups were found. The OS rates were also compared, and we found significant differences (p = 0.013, Figure 3B ). To study the prognostic effects of these genes, the PCA algorithm was also applied to determine the PCA gene score of each individual, and cases were divided into two groups. The KM curves showed that the high-score group had obviously a better prognosis (p < 0.001, Figure 3C ).

Based on the above analysis, this study used ssGSEA to compare the enrichment fraction of 16 types of immune cells and the activity difference of 13 immune-related pathways among groups in detail. The comparative analysis results show that in the TCGA queue ( Figures S2A, B ), the immune cell infiltration of the subgroup with high gene score is also higher, and the corresponding immune infiltration is mainly related to activated CD8⁺T, CD4⁺T, and Treg cells. The comparison results also showed that the infiltration of neutrophils and NK and Th cells in the high-score subgroup was more obvious.

Some immune pathways, such as the TGF beta pathway, NK cell-mediated cytotoxicity pathway, JAK–STAT pathway, Toll-like receptor pathway, and B-cell receptor pathway, showed higher activity in the high-score group ( Figure S2C ). Moreover, amino sugar/nucleotide sugar/starch/sucrose metabolism and other glycan degradation pathways have higher activity in the high-gene score group, while taurine/hypotaurine-glycine/serine/threonine metabolism, calcium signaling, and neuroactive ligand–receptor interactions have lower activity in the high-gene score group ( Figure S2D ). According to the comparative analysis results, the activity of the sucrose degradation pathway in the high-gene score group is significantly higher than that in the other group, but the corresponding calcium signaling and neural activity are relatively low ( Figure S2D ). The cause is not very clear, and it needs to be analyzed in detail to provide support for clarifying the pathological correlation between ferroptosis and this disease, so as to better meet the requirements related to treatment. Correlation analysis showed a significant positive correlation between immune cells in GBM, and almost all immune cells have a positive association with the PCA score and PCA gene score ( Figures 4A, B ).

In order to further verify the role of the 25 differential genes we screened earlier in glioma, we detected the expression of the above genes in glioma glial tissue and glioma paracancerous glial tissue by immunohistochemistry, and the expression of the above genes in glioma glial tissue was significantly increased. Next, we detected four genes with differential expression greater than 10-fold by Western blot at the cellular and tissue levels, namely, MAP1LC3A, OLFM1, CEND1, and CLTB. The results showed that the expression levels of the above proteins were significantly upregulated relative to both glioma paracancerous glial cells and glioma paracancerous glial tissue. This further validates the results of our previous analysis ( Figure 5 ).