MRNA expression data, DNA methylation 450K data, and corresponding clinical information of 529 LGG patients were obtained from TCGA data portal (https://gdc-portal.nci.nih.gov/). For differentially expressed gene analysis, gene expression data in normal tissue of Genotype-Tissue Expression (GTEx) were downloaded from the UCSC Xena browser (Goldman et al., 2014; Gentles et al., 2015). RNA-seq data and clinical information of other LGG sets were obtained from the Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn/) (Bao et al., 2014; Wang et al., 2015; Zhao et al., 2017; Liu X. et al., 2018). Gene expression values of fragments per kilobase of exon per million reads mapped (FPKM) were normalized and subjected to the following analyses. A comprehensive list of ferroptosis-related genes (FRGs) was chosen and is provided in Supplementary Table S1.

To identify differentially expressed FRGs (DE-FRGs) between TCGA and GTEx databases, the absolute value of the log₂-transformed fold change (log₂FC) > 2 and the false discovery rate (FDR) < 0.05 was considered to be statistically significant (Zhang et al., 2018). To investigate the mechanisms of dysregulation for FRGs in LGG, Pearson’s correlation between FRGs’ transcriptional expression levels and promotor methylation levels was calculated, and those with |R| > 0.3 and p < 0.05 were considered significantly correlated (Xu et al., 2019).

To screen genes for constructing a risk signature, univariate Cox regression models (p < 0.05) were performed to select genes that are associated with the overall survival of LGG patients in TCGA. The overlapped candidate prognostic DE-FRGs were visualized in a Venn diagram. A protein–protein interaction (PPI) network for the prognostic DE-FRGs was generated by the STRING database (version 11.0) (Szklarczyk et al., 2010).

Subsequently, the least absolute shrinkage and selection operator (LASSO) Cox regression was performed by using the “glmnet” R package to shrink the scope of gene screening. The multivariate Cox regression analysis was carried out to identify highly correlated genes and construct the survival risk score model. The risk score was calculated as follows: Risk   score= ∑ i = 1 n expi × βi , where n, expi, and βi represent the number of genes, the expression level, and the coefficient of the gene i, respectively. According to the median value of the risk scores, the LGG patients were classified into low- and high-risk groups. Principal component analysis (PCA) was carried out to explore the distribution of different groups using the “prcomp” function in R. The Kaplan–Meier (K–M) curves were used to evaluate the statistical significance of the survival rates between different risk groups. The receiver-operating characteristic (ROC) curve analyses were executed to calculate the area under the curve (AUC) and to determine the predictive ability of the risk score signature. We calculated the false negative rate (FNR) of the survival prognosis by our FRG model.

In order to study whether the FRG-based prognostic signature could be used as an independent predictor for LGG patients, univariate and multivariate Cox regression analyses were conducted in TCGA and CGGA datasets. The risk score, age, gender, grade, and IDH status were used as covariates.

To establish a clinically applicable method for predicting the prognosis of LGG patients, we formulated a prognostic nomogram to predict the survival probability at 1, 3, and 5 years based on TCGA dataset. Five independent prognostic parameters, including age, gender, grade, IDH1 status, and risk score, were used to establish the nomogram for predicting one-, three-, and five-year survival rates of LGG patients in the “rms,” “foreign,” and “survival” R packages. Bootstraps with 1,000 resamples were used for these activities. The performance of the prognostic nomogram was assessed by calculating Harrell’s concordance index (C-index) (Harrell et al., 1996), and calibration curves of the nomogram for one-, three-, and five-year overall survival were carried out to estimate the accuracy of actual observed rates with the predicted survival probability. During the external validation of the nomogram, the total points of each patient in the CGGA cohort were calculated according to the established nomogram, and then the C-index was calculated and calibration curves were derived.

The “clusterProfiler,” “org.Hs.eg.db,” “enrichplot,” and “ggplot2” packages in R were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation enrichment analyses based on the differentially expressed genes (|log₂FC| ≥ 2, FDR <0.05) between the low- and high-risk groups. The top 30 GO terms and KEGG pathways were identified with a cutoff of p < 0.05 and FDR < 0.05.

The single-sample gene set enrichment analysis (ssGSEA) was conducted using the “GSVA” package in R to evaluate the enrichment score (ES) of immune parameters, including 15 immune cells and 13 immune-related pathways (Goldman et al., 2014). The annotated gene set file is provided in Supplementary Table S2. Pearson’s correlation coefficients were calculated to test the correlation between the ES and the risk score.

All statistical analyses were performed using R (version 4.0.3). Differences between two sample groups were tested by the Wilcoxon test. Univariate and multivariate Cox regression analyses were performed to determine risk factors for low and high risk scores of the signature in LGG. The Venn diagrams, heatmaps, boxplots, forest plots, nomograms, and calibration plots were drawn using R language. The K–M curve used the log-rank test to evaluate the statistical difference of survival rates between different risk groups. The significance of survival time differences was calculated using the log-rank test with a threshold of p < 0.05.

The workflow of this study is shown in Figure 1. A total of 508 LGG patients from TCGA-LGG cohort and 171 LGG patients from the CGGA cohort were finally enrolled. The clinical characteristics of patients from two databases are summarized in Table 1.

32 FRG expressions were identified to have a significant difference between LGG samples from TCGA and normal samples from GTEx (Figure 2A). Fifteen genes were upregulated in LGG, whereas another seventeen genes were downregulated in tumor. Since methylation of promotor plays a critical role in regulating gene expression, the effect of promotor DNA methylation on the gene expression of DE-FRGs was evaluated. As shown in Figures 2A,B, a significant negative correlation was found between gene expression and DNA methylation in five genes, including FTH1, ACSL3, SLC1A5, HSPB1, and STEAP3, whereas significant positive correlations were found in CRYAB and ZEB1. These results demonstrated that promoter DNA methylation might regulate the expression of FRGs in LGG with complexly regulative patterns.

A total of 37 FRGs were found for potent prognostic value (p < 0.05) through univariate Cox proportional hazard regression. 20 of them were differentially expressed between LGG and normal tissue (Figures 2C,D). The PPI network among these genes was constructed using the STRING online database, with parameters including a minimum required interaction score >0.4 (medium confidence). As shown in Figure 2E, TP53, SLC7A11, and HMOX1 were hub genes for ferroptosis in LGG (Figure 2E). LASSO regression analysis was performed on the 20 overlapped genes, and nine genes were retained according to the optimal lambda value (Supplementary Figure S1). By using multivariate Cox regression analysis, a prognostic model based on seven genes was developed. The prognostic risk score was imputed as follows: Risk score = (0.0086 × expression level of ACSL3) + (−0.0190 × expression level of CBS) + (0.0013 × expression level of CD44) + (−0.0027 × expression level of FADS2) + (0.0017 × expression level of HSPB1) + (0.0104 × expression level of PGD) + (0.0271 × expression level of STEAP3). The 509 patients were classified into high- and low-risk groups according to the median value as cutoff (Figure 3A and Supplementary Figures S2A,C). The PCA plot showed the patients in two risk groups were distributed in discrete directions (Figure 3B). Kaplan–Meier (K–M) plots indicated patients with high risk scores had a significantly worse OS probability (p = 9.178e-9) (Figure 3C). To verify the diagnostic competence of the ferroptosis-related risk signature, the AUC for OS prediction was calculated. The AUC of the ROC curve reached 0.855 at 1 year, 0.834 at 3 years, and 0.753 at 5 years, which indicated that the risk score presented excellent performance on prognostic prediction of TCGA-LGG (Figure 3D). The false negative rate (FNR) of prognosis for OS at 1 , 3 , and 5 years was 0.077, 0.118, and 0.169, respectively (Supplementary Table S3).

The prognostic capability of the FRG risk signature was externally validated by the CGGA dataset. The patients from the CGGA cohort were also divided into low- or high-risk groups by the same cutoff value as that in TCGA cohort (Figure 4A, Supplementary Figures S2B,D, and Supplementary Figure S3). Consistent with the above findings, PCA found that the patients in different risk groups were distributed in two directions (Figure 4B). Additionally, the K–M survival curves showed a significantly worse survival rate with high-risk patients (Figure 4C and Supplementary Figure S3A). The AUC for one-, three-, and five-year survival was 0.752, 0.795, and 0.794, respectively (Figure 4D). These results confirmed that the seven-FRG prognostic model presented good performance in prognostic prediction of LGG.

To evaluate the independent prognostic force of the signature, both the univariate and multivariate Cox proportion hazard regression models were applied. The risk score in univariate analysis was significantly correlated with OS in both TCGA and CGGA cohorts (HR = 3.814, 95% CI = 2.021–7.196, p < 0.001; HR = 4.250, 95% CI = 2.719–6.644, p < 0.001, respectively) (Figures 5A,B). The results from multivariate analysis found this signature was an independent prognostic factor for OS in TCGA-LGG (HR = 2.352, 95% CI = 1.179–4.690, p = 0.015) and CGGA-LGG (HR = 2.758, 95% CI = 1.572–4.840, p < 0.001) cohorts. These results confirmed the independent prognostic value of the seven-FRG–based prognostic signature for patients with LGG (Figures 5C,D).

Next, by integrating the above risk signature with clinicopathological factors, including age, gender, grade, and IDH status, a nomogram was cautiously constructed for prediction of one-, three-, and five-year survival rates of LGG. As shown in Figure 6A, each level of factors was assigned one point, and a total point was calculated by adding up the points in each individual. Subsequently, the C-indexes of the nomogram for predicting OS were 0.87 (95% CI, 0.83–0.91) and 0.68 (95% CI, 0.63–0.73) in TCGA and CGGA cohorts, respectively. Calibration curves also showed good agreement between predictive and observational values at the probabilities of one-, three-, and five-year survival in TCGA and CGGA cohorts (Figures 6B,C). These results revealed that the nomogram signified good accuracy in predicting the one-, three-, and five-year survival of patients with LGG.

To further elucidate the biological functions and pathways associated with the risk signature, we performed GO and KEGG pathway enrichment analyses based on differentially expressed genes (DEGs) between the high- and low-risk groups. As shown in Figures 7A,C, GO analysis showed 13 clusters sets were enriched in both TCGA and CGGA cohorts. For biological processes (BPs), DEGs were mainly enriched into immune-related terms, including neutrophil degranulation and neutrophil activation involved in immune response. For cellular components (CCs), cell-substrate junction, focal adhesion, cell leading edge, and secretory granule lumen were found to be enriched. With respect to molecular function (MF), differentially expressed genes were mainly enriched in actin binding, integrin binding, cadherin binding, actin filament binding, and growth factor binding. For the KEGG pathway, 10 out of top 30 pathways were enriched in both cohorts, and most of these enriched pathways were related to infections of bacterial or viral (Figures 7B,D). Interestingly, the term of focal adhesion, which was also reported closely related to the tumor immune microenvironment (TME) (Murphy et al., 2020), was enriched in both GO and KEGG analyses. These results revealed that the risk signature was correlated with the infection and immune-related pathway.

Furthermore, we quantified the enrichment scores (ESs) of diverse immune cell subtypes and related functions or pathways with ssGSEA to confirm the correlation between the risk signature and immune status. As shown by the boxplot in Figure 8A, in TCGA-LGG cohort, 25 out of 28 immune-related terms were significantly higher in the high-risk group (p < 0.05, Figure 8A), except DCs (dendritic cells), neutrophils, and NK cells (natural killer cells). Consistently, twenty-four immune-related terms were validated by the ICGC cohort, except Tfh cells (T follicular helper cells), while NK cells were significantly higher in the high-risk group (p < 0.01, Figure 8C).

We further investigated the association between expressions of the seven-gene signature and the immune status ES in detail by calculating Pearson’s correlation. The results showed that most of the seven genes were positively correlated with these immune cells and functions, but CBS (cystathionine beta-synthase) and FADS2 (fatty acid desaturase 2) were negatively correlated with these immune terms in both cohorts (p < 0.05, Figures 8B,D). However, inconsistent correlations of acyl-CoA synthetase long-chain family member 3 (ACSL3) with immune characters were observed between TCGA and CGGA cohorts (Figures 8B,D).