Iron metabolism-related genes were retrieved from gene sets downloaded from the Molecular Signatures Database (MSigDB) version 7.1 (16, 17), including the GO_IRON_ION_BINDING, GO_2_IRON_2_SULFUR_CLUSTER_BINDING, GO_4_IRON_4_SULFUR_CLUSTER_BINDING, GO_IRON_ION_IMPORT, GO_IRON_ION_TRANSPORT, GO_IRON_COORDINATION_ENTITY_TRANSPORT, GO_RESPONSE_TO_IRON_ION, MODULE_540, GO_IRON_ION_HOMEOSTASIS, GO_CELLULAR_IRON_ION_HOMEOSTASIS, GO_HEME_BIOSYNTHETIC_PROCESS, HEME_BIOSYNTHETIC_PROCESS, GO_HEME_METABOLIC_PROCESS, HEME_METABOLIC_PROCESS, HALLMARK_HEME_METABOLISM, and REACTOME_IRON_UPTAKE_AND_TRANSPORT gene sets. We also reviewed the literature and added the previously reported genes (18, 19). After removing overlapping genes, we obtained an iron metabolism-related gene set containing 527 genes.

Gene expression data for 523 LGG samples (TCGA) and 105 normal cerebral cortex samples (GTEx project) were downloaded from a combined set of TCGA, TARGET, and GTEx samples in UCSC Xena (https://tcga.xenahubs.net). Clinical information for patients with LGG was obtained from using the “TCGAbiolinks” package written for R (20–22). Gene expression data and clinicopathological information for 443 patients with LGG were retrieved from CGGA database (http://www.cgga.org.cn/) and were selected as a test set. Data from patients without prognostic information were excluded from our analysis. Ultimately, we obtained a TCGA training set containing 506 patients and a CGGA test set with 420 patients. Ethics committee approval was not required since all the data were available in open-access format.

First, we screened out 402 duplicate iron metabolism-related genes that were identified in both TCGA and CGGA gene expression matrixes. Then, differentially expressed genes (DEGs) between the TCGA-LGG samples and normal cerebral cortex samples were analyzed using the “DESeq2”, “edgeR” and “limma” packages of R software (version 3.6.3) (23–26). The DEGs were filtered using a threshold of adjusted P-values of < 0.05 and an absolute log₂-fold change >1. Venn analysis was used to select overlapping DEGs among the three algorithms mentioned above. Eighty-seven iron metabolism-related genes were chosen for downstream analyses. Additionally, functional enrichment analysis of selected DEGs was performed using Metascape (https://metascape.org/gp/index.html#/main/step1) (27).

Univariate Cox proportional hazards regression was performed for the genes selected for the training set using “ezcox” package (28). P < 0.05 was considered to reflect a statistically significant difference. To reduce the overfitting high-dimensional prognostic genes, the Least Absolute Shrinkage and Selection Operator (LASSO)-regression model was performed using the “glmnet” package (29). The expression of identified genes at protein level was studied using the Human Protein Atlas (http://proteinatlas.org). Subsequently, the identified genes were integrated into a risk signature, and a risk-score system was established according to the following formula, based on the normalized gene expression values and their coefficients. The normalized gene expression levels were calculated by TMM algorithm by “edgeR” package.

The risk score was calculated for each patients with LGG in this study, and the distribution and receiver operating characteristic (ROC) curve were plotted using “timeROC” package (30). According to the median risk score in the training set, patients were divided into high- or low- risk groups. Patients were also divided into subgroups according to clinicopathological features, including age, gender, WHO grade, histological type, isocitrate dehydrogenase 1 (IDH1) mutation status, 1p19q codeletion status, and O⁶‐methylguanine‐DNA methyl‐transferase (MGMT) promoter methylation status. Boxplot were plotted using the “ggpubr” package to identify associations between risk scores and clinical features. In addition, the relationships between risk scores and OS rates in different groups and subgroups were evaluated by Kaplan-Meier survival analysis and log-rank testing.

To evaluate whether the risk score system can serve as an independent predictive index, univariate and multivariate Cox regression analyses were performed with clinicopathological parameters, including the age, gender, WHO grade, IDH1 mutation status, 1p19q codeletion status, and MGMT promoter methylation status. All independent prognostic parameters were used to construct a nomogram to predict the 1-, 3- and 5-year OS probabilities by the ‘rms’ package. Concordance index (C-index), calibration and ROC analyses were used to evaluate the discriminative ability of the nomogram (31).

DEGs between high- and low-risk groups in the training set were calculated using the R packages mentioned above. Then, GSEA (http://software.broadinstitute.org/gsea/index.jsp) was performed to identify hallmarks of the high-risk group compared with the low-risk group.

The TIMER database (http://timer.cistrome.org/) is a comprehensive web tool that provide automatic analysis and visualization of immune cell infiltration of all TCGA tumors (32, 33). The infiltration estimation results generated by the TIMER algorithm consist of 6 specific immune cell subsets, including B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils and dendritic cells. We extracted the infiltration estimation results and assessed the different immune cell subsets between high-risk and low-risk groups (34).

All statistical analyses in this study were conducted using R software (version 3.6.3) and GraphPad Prism (version 8.0.2). The log-rank test was used for the Kaplan-Meier survival analysis. Hazard ratios (HRs) and 95% confidence intervals (CIs) were reported where applicable. Student’s t-test and the Kruskal–Wallis test were employed in the two-group comparisons. A two-tailed P value of <0.05 was considered statistically significant without specific annotation.

The data we used were retrieved from open-access databases. The majority of statistical codes are available in File S1.

Based on the MSigDB and a literature review, we selected 527 iron metabolism-related genes for analysis. Four hundred and two genes remained after excluding genes not present in the TCGA-LGG or CGGA-LGG set. According to the criteria for DEG, we identified 7,223 DEGs between 523 TCGA-LGG samples and 105 normal brain cortex samples based on overlapping edgeR, limma, DESeq2 analysis results ( Figure 1A ). Then, a total of 87 iron metabolism-related genes (50 up-regulated and 37 down-regulated) among the DEGs were selected for further analysis ( Figure 1B ). Enrichment analyses were performed to explore the functions of the selected genes. These genes were significantly enriched in terms of iron ion binding, iron ion metastasis, and iron ion transport ( Figure 1C ). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that ferroptosis, mineral absorption, the p53 signaling pathway and the AMPK signaling pathway were enriched ( Figure 1D ).

First, univariate Cox regression was used to investigate the relationship between the expression levels of the selected genes and OS time in the training set. Using cut-off threshold of Cox P < 0.05, 47 genes were identified as potential risk factors related to OS ( Table S1 ). Subsequently, the LASSO regression algorithm was used to refine the gene sets by calculating regression coefficients ( Figures 2A, B ). In this manner, 15 genes were identified as the most valuable predictive genes, and the risk-score system was established using the formula mentioned above ( Table 1 ).

We also confirmed the expression level of these identified genes by Immunohistochemical analysis in Human Protein Atlas (HPA). And the results were shown in Figure 3 . 6 of these genes were dysregulated in LGG and higher-grade glioma samples. The expressions level of GCLC, NCOA4, UROS were higher in LGG samples, whereas the expression levels of LAMP2, RRM2, STEAP3 were lower in LGG than HGG samples. CH25H and RTEL1 were missing in HPA database. ACP5, CYP2D6, HBQ1, KHNYN, and SCD5 were not detected in glioma samples. However, the expression levels of CYP2E1 and FLVCR2 showed low consistency with RNA expression data.

The risk score for each patient in the training and test sets was calculated based on the expression levels of the selected genes and the regression coefficients. The distribution of risk score in training set was shown in Figure 4A . The median of risk score in training set was defined as threshold, which divided the patients into high-risk and low-risk groups. In addition, the distribution of survival times indicated that a higher risk score may have positively correlated with poorer outcomes ( Figure 4A ). The corresponding expression levels of the selected genes were determined ( Figure 4A ). The performance of the ROC in terms of 1-, 3-, and 5-year prognoses was analyzed ( Figure 4B ). The areas under the time‐dependent ROC curve (AUCs) were 0.892, 0.888, and 0.838, respectively, for the 1-, 3-, and 5-year OS times in the training set. Kaplan–Meier analysis and log-rank testing showed that the high-risk group had a significantly shorter OS time than the low-risk group (P < 0.0001; Figure 4C ).

Furthermore, the robustness of our risk-score model was assessed with the CGGA dataset. The test set was also divided into high-risk and low-risk groups according to the threshold calculated with the training set. The distributions of risk scores, survival times, and gene-expression level are shown in Figure 4D . The AUCs for the 1-, 3-, and 5-year prognoses were 0.765, 0.779, and 0.749, respectively ( Figure 4E ). Significant differences between two groups were determined via Kaplan–Meier analysis (P < 0.0001), indicating that patients in the high-risk group had a worse OS ( Figure 4F ). These results showed that our risk score system for determining the prognosis of patients with LGG was robust.

Associations between risk-score and clinical features in the training set were examined. We found that the risk score was significantly lower in groups of patients with age > 40 (P < 0.0001), WHO II LGG (P < 0.0001), oligodendrocytoma (P < 0.0001), IDH1 mutations (P < 0.0001), MGMT promoter hypermethylation (P < 0.0001), and 1p/19q co-deletion (P < 0.0001) ( Figures 5A–F ). However, no difference was found in the risk scores between males and females (data not shown). In both astrocytoma and oligodendrocytoma group, risk score was significantly lower in WHO II group ( Figures 5G, H ). We also validate the prediction efficiency with different subgroups. Kaplan–Meier analysis showed that high-risk patients in all subgroups had a worse OS ( Figure S1 ). Besides, the risk score was significantly higher in GBM group compared with LGG group ( Figure S2 ).

To determine whether the risk score was an independent risk factor for OS in patients with LGG, the potential predictors (age group, gender, WHO grade, IDH1 mutation status, MGMT promoter status, 1p/19q status and risk level) were analyzed by univariate Cox regression with the training set ( Table 2 ). The individual risk factors associated with a Cox P value of < 0.05 were further analyzed by multivariate Cox regression ( Table 2 ). The analysis indicated that the high-risk group had significantly lower OS (HR = 2.656, 95% CI = 1.51-4.491, P = 0.000268). The age group, WHO grade, IDH mutant status, MGMT promoter status and risk level were considered as independent risk factors for OS, and were integrated into the nomogram model ( Figure 6A ). The C-index of the nomogram model was 0.833 (95% CI = 0.800-0.867). Subsequently, we calculated the score of each patient according to the nomogram, and the prediction ability and agreement of the nomogram was evaluated by ROC analysis and a calibration curve. In the TCGA cohort, the AUCs of the nomograms in terms of 1-, 3-, and 5-year OS rates were 0.875, 0.892, and 0.835, respectively ( Figure 6B ). The calibration plots showed excellent agreement between the 1-, 3-, and 5-year OS rates, when comparing the nomogram model and the ideal model ( Figures 6D–F ). Moreover, we validated the efficiency of our nomogram model with the CGGA test set. The AUCs for the 1-, 3-, and 5-year OS rates with the model were 0.722, 0.746, 0.701, respectively ( Figure 6C ). The results of the calibration curves showed good agreement between the predicted OS rates and the probabilities of the 1-, 3-, and 5-year OS rates with the test set ( Figures 6G–I ).

To clarify the potential impact of the expression levels of the selected genes on the LGG transcriptomic profile, GSEA analysis was performed with the high-risk and low-risk groups of the training set. GSEA revealed that several pathways, such as those related to inflammatory response, IL6/JAK/STAT3 signaling, IL2/STAT5 signaling, glycolysis, apoptosis, and the EMT, were enriched in the high-risk group ( Figures 7A–F ). These findings suggest potential roles for iron metabolism-related genes in the progression, metabolism, tumor microenvironment and immune responses of LGG.

Next, the correlation between this prognostic model and the infiltration of immune cells for patients in the TCGA-LGG cohort were calculated. The proportion of different infiltrating immune cells were retrieved from the TIMER database. The results indicated that the risk score positively correlated with infiltrating immune cells, including B cells, CD4+ T cells, CD8+T cells, neutrophils, macrophages and dendritic cells ( Figure 8A ). The high-risk group showed more infiltrating immune cells, especially dendritic cells and macrophages (P < 0.0001; Figure 8B ). Additionally, we assessed the relationship between risk-score model and immune checkpoint proteins (PD1, PDL1, CTLA4, LAG-3, TIM3, TIGIT and CD48). The expression levels of PD1, PDL1, CTLA4, TIM3, and CD48 positively correlated with the risk score(P < 0.001; Figure 8C ). In addition, the expression levels of PD1, PDL1, and TIM3 were higher in high-risk group of TCGA-LGG cohort than in the low-risk group (P < 0.0001; Figure 8D ).