Ferroptosis‐related gene signature as a prognostic marker for lower‐grade gliomas

Abstract Ferroptosis is a newly discovered form of programmed cell death, which has unique biological effects on metabolism and redox biology. In this study, the prognostic value of ferroptosis‐related genes was investigated in lower‐grade gliomas (LGG). We downloaded the ferroptosis‐related genes from the FerrDb dataset. Univariate Cox and LASSO regression analyses were applied to identify genes correlated with overall survival (OS). Subsequently, 12 ferroptosis‐related genes were screened to establish the prognostic signature using stepwise multivariate Cox regression. According to the median value of risk scores, patients were divided into low‐ and high‐risk subgroups. The Kaplan‐Meier curves showed the high‐risk group had a lower OS. The predictive power of the risk model was validated using the CGGA. Functional analysis revealed that the terms associated with plasma membrane receptor complex, immune response and glutamate metabolic process were primarily related to the risk model. Moreover, we established a nomogram that had a strong forecasting ability for the 1‐, 3‐ and 5‐year OS. In addition, we compared the risk scores between different clinical features. We also detected infiltration of macrophages and monocytes in different subgroups. Overall, our study identified the prognostic signature of 12 ferroptosis‐related genes, which has the potential to predict the prognosis of LGG.

associated with the imbalance of redox homeostasis. 7 The cystine/glutamate antiporter (system Xc-) is a redox state regulator on the surface of the cellular plasma membrane, 8 which can regulate the exchange of intracellular glutamate and extracellular cystine. 9 Cystine promotes the synthesis of glutathione (GSH), which plays an important role as an intracellular antioxidant. 10 Certain molecules deplete GSH by inhibiting the function of system Xc-, leading to the accumulation of irondependent lipid-ROS. 11 Tumour cells, enriched in free iron and with high levels of ROS, are more sensitive to ferroptosis. Recent studies have reported the association of ferroptosis with glioma cell proliferation, invasion and angiogenesis. [12][13][14] However, there are few studies on its correlation with the prognostic value in patients with LGG.
In this study, we collected the RNA-Seq data and clinical information of LGG in The Cancer Genome Atlas (TCGA). We identified 12 ferroptosis-related genes by statistical analysis to establish a prognostic risk model. Meanwhile, patients with LGG in the Chinese Glioma Genome Atlas (CGGA) database were selected as a validation cohort. Gene Ontology (GO), 15,16 Kyoto Encyclopedia of Genes and Genomes (KEGG), 17 Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA) were used to screen the functions and pathways enriched between high-risk and low-risk subgroups.
Furthermore, we developed a nomogram model based on risk scores and clinical features to assess prognosis. Our results demonstrated that the ferroptosis-related risk model was a potential prognostic marker and therapeutic target for LGG. The overview workflow is presented in Figure 1.

| Download of datasets and collection of ferroptosis-related genes
The gene expression RNA-Seq (HTSeq-FPKM) and clinical information of 495 patients with LGG in TCGA were collected from the University of California Santa Cruz (UCSC) Xena (https://xenab rowser.net/datap ages/). Furthermore, 590 LGG samples were downloaded from the CGGA (http://cgga.org.cn/). In addition, the exclusion criteria included patients with a follow-up of less than 30 days and a single gene with a total expression value of less than 10 in all 495 samples. We downloaded 259 ferroptosis-related genes from the FerrDb data set (http://www.zhoun an.org/ferrd b/).

| Construction and validation of the prognostic risk model
A univariate Cox regression analysis was used to screen genes associated with prognosis, and we considered P < .01 as statistically significant. The least absolute shrinkage and selection operator (LASSO) regression was used to reduce the potential risk for overfitting with the 'glmnet' package in R. 18 Ten-fold cross-validation (10FCV) was performed to select the optimal value of λ, and 22 key genes were obtained. We then used the stepwise multivariate Cox F I G U R E 1 Flowchart of the construction of the ferroptosis-related genes prognostic signature regression to perform the prognostic signature with 12 prognostic genes and drew a forest map to visualize the result.

The formula
where n represents the total number of genes, β i is the coefficient, and X i is the expression value of the gene, was used to calculate the risk score: risk score = (β 1 *expression of gene X 1 ) + (β 2 *expression of gene X 2 ) + … (β 12 *expression of gene X 12 ). Patients were clustered into low risk and high risk based on the median value of risk score in TCGA. Next, CGGA was used to verify the predictive power of the risk model. The Kaplan-Meier survival curves and log-rank tests were used to compare the differences in overall survival (OS) between the two categories. The risk and survival status plots are shown by the rank of the corresponding risk scores ( Figure 2E,H). Finally, the 1-, 2-, 3-and 5-year receiver operating characteristic (ROC) curves were plotted for TCGA and CGGA datasets to evaluate the sensitivity and specificity of survival prediction using 'timeROC' in R.

| Identification and functional enrichment of differentially expressed genes
The expression of 12 signature genes and their corresponding clinical information are shown in the heat maps. We identified the differentially expressed genes (DEGs) using the Wilcoxon test between low-and high-risk groups with FDR < 0.05 and |logFC| > 1. DEGs were analysed using the 'clusterProfiler' package in R for GO and KEGG. GSEA 19 was conducted to confirm the activation or inhibition of biological processes and signalling pathways in low-and high-risk groups. The KEGG gene set (C2.cp.kegg.v7.2.symbols.gmt) and GO gene set (C5.go.v7.2.symbols.gmt) were adopted from the Molecular Signatures Database (MSigDB). An FDR q < 0.25 and adjusted P < .05 were considered statistically significant. As a non-parametric unsupervised analysis method, GSVA 20 was used to assess the gene set enrichment (GSE) for each sample. We applied the 'GSVA' package in R to investigate significant metabolic and immunologic pathway differences between groups in the TCGA and CGGA datasets.

| Nomogram construction and evaluation
Univariate and multivariate Cox analyses were used to screen for prognostic factors such as grade, radiation therapy, age, gender, IDH1 status, MGMT promoter status, 1p/19q codeletion and risk score. We filtered variables by stepwise regression to avoid the effect of multicollinearity. 21 Subsequently, multivariate Cox regression analysis was used to identify the independent prognostic factors of the model. Thereafter, a nomogram was constructed, using the result of multi-Cox analysis, to predict the 1-, 3-and 5-year OS in patients with LGG. This was then assessed by calibration curves. In addition, we compared the risk scores between different clinical features.

| Estimation of the abundance of macrophages and monocytes by ImmuCellAI
ImmuCellAI (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) can estimate the abundance of immune cell infiltration with higher accuracy than other algorithms. 22 Using the ImmuCellAI website, we contrasted the infiltration levels of macrophages and monocytes between low-and high-risk subgroups in TCGA and CGGA. P < .05 were accepted as statistically significant.

| Construction of the ferroptosis-related prognostic risk model
We searched for genes associated with ferroptosis in the FerrDb dataset. A total of 241 ferroptosis-related genes were extracted in the TCGA-LGG cohort. Univariate Cox and LASSO regression analyses were applied to screen for 22 feature genes (Figure 2A,B).
Then, stepwise multivariate Cox regression analysis was used to select the best characteristic gene set and construct a regression model. Finally, the prognostic signature was established based on 12 ferroptosis-related genes ( Figure 2C). Their functions and coefficients are shown in Table S1. The cut-off value of risk scores (−0.344) was used to dichotomize patients into low-risk (n = 248) and high-risk (n = 247) groups in the TCGA cohort. The prognosis was significantly better in the low-risk group than in the high-risk group (P <.001, HR = 0.19, 95% CI 0.12-0.28; Figure 2D). As shown in Figure 2E, the distributions of risk score and survival status indicated that low-risk scores were advantageous for survival. The ROC curves confirmed that the model had a good accuracy for predicting OS in TCGA (1-year AUC = 0.902, 2-year AUC = 0.919, 3-year AUC = 0.925, 5-year AUC = 0.837; Figure 2F).

| Validation of the prognostic signature in the CGGA cohort
We examined the predictive power of the prognostic signature in the CGGA data set. It was divided into low-risk (n = 297) and highrisk (n = 293) groups using the same formula and cut-off value as those for TCGA. As shown in Figure 2G

| Functional analysis of risk model in TCGA and CGGA
The heat map of 12 genes in the ferroptosis-related prognostic risk model, combined with clinical information, illustrated that patients in the high-risk group had a tendency towards higher expression levels of risk genes and lower expression levels of protective genes ( Figure 3A). These results were similar to those for CGGA ( Figure 3B). A total of 1111 DEGs were identified between the highand low-risk groups in the TCGA cohort. GO analysis was primarily associated with the plasma membrane receptor complex, cytokine and immune response-related terms ( Figure 3C). KEGG analysis revealed that cytokine-cytokine receptor interaction, T cell receptor signalling pathway, Th1 and Th2 cell differentiation, and tyrosine metabolism were obtained ( Figure 3E). Furthermore, 806 DEGs were selected for GO and KEGG in the CGGA cohort. The results revealed that immune response, cytokine activity, receptor ligand activity and neuroactive ligand-receptor interaction were primarily enriched ( Figure 3D,F).
GSEA analysis showed that the high-risk group was associated with the pathways of glutathione metabolism and immune-related functions from the KEGG database ( Figure 4A). Among GO terms, the glutamate metabolic process, glutamate receptor binding, glutamate receptor activity, glutamate receptor signalling pathway and ionotropic glutamate receptor signalling pathway, which play vital roles in ferroptosis, were found to be enriched in the low-risk group ( Figure 4B). GSVA was used to identify the causes of the phenotypic differences. We found that the glutathione disulfide oxidoreductase activity, glutathione peroxidase activity, ferric iron binding, cellular response to iron ion, negative regulation of T cell differentiation, regulation of iron ion transmembrane transport and negative regulation of T cell-mediated cytotoxicity were enriched in the high-risk group. Furthermore, type 5 metabotropic glutamate receptor binding, glutamate receptor activity and AMPA glutamate receptor clustering were more enriched in the low-risk group ( Figure 4E). These results were validated in CGGA ( Figure 4C,D,F).

| Development and validation of a nomogram
Univariate and multivariate Cox regression analyses were performed to identify the risk score as an independent prognostic biomarker in LGG ( Figure 5A,B). The nomogram used five prognostic markers  Figure 5E). The calibration plot of validation set showed that the predicted power was close to the ideal curve ( Figure 5F). In addition, the C-index of the nomogram was 0.817 (95% CI 0.762-0.872). Overall, these results demonstrated that the developed nomogram preforms well in predicting OS.

| Correlation between clinical features and risk score
There were significant correlations between clinical parameters (age, grade, IDH1 status, MGMT promoter status, radiation therapy and 1p/19q codeletion) and the risk score. We found that the risk score was higher in older patients (age > 40 years), and in patients with grade III, IDH1 wild-type, MGMT promoter unmethylated, patients receiving radiation therapy, and in those with 1p/19q noncodeletion (P <.05, Figure 6A-F).

| Immune cell infiltration in high-and low-risk groups
We assessed the differences in macrophage and monocyte infiltration between different subgroups, as defined by the risk model. Figure 6G shows that the high-risk group had significantly higher populations of macrophages and monocytes than the low-risk group.
As for the CGGA database, the results were similar to those for the TCGA ( Figure 6H). It is worth mentioning that ferroptosis may regulate the tumour immune microenvironment, which has great significance for glioma therapy.

| D ISCUSS I ON
Increasing evidence has shown that ferroptosis is essential for eradicating carcinogenic cells 23 and that the sensitivity to ferroptosis is different in various types of cancers. 24 Cysteine production, lipid peroxidation of polyunsaturated fatty acids (PUFAs), iron metabolism and mitochondrial function are closely related to ferroptosis.
Activation of Nrf2-Keap1 signalling up-regulates the glutamatecystine antiporter system, which accelerates the progression of glioma. 25 However, the prognostic value of ferroptosis-related genes in LGG has yet to be clarified. In this study, we identified 12 ferroptosisrelated genes to construct a prognostic risk model. Furthermore, the risk score was an independent prognostic biomarker in patients with LGG, which was closely related to the clinical features and immune functions.  Jiang et al found that DNAJB6 plays an important role in ferroptosis, and its expression is down-regulated in the oesophageal carcinoma tissues. 27 In addition, the overexpression of DNAJB6 leads to radiosensitization of glioblastoma cells. 28 Among the other eight risk-associated genes, NFE2L2, also known as Nrf2, mediated oxidative stress and inflammatory response. 29,30 Activation of Nrf2-Keap1 It has a cytoprotective effect on glioma through the inhibition of apoptosis. 32 HSPB1 correlated with poor outcomes and promoted the proliferation of glioma cells by facilitating an anti-oxidative response. 33 Previous studies have shown that HSPB1 is a negative regulator of the erastin-mediated ferroptosis. 34 ARNTL/BMAL1 is associated with molecular circadian rhythms. A recent study demonstrated that the inhibition of the ARNTL-EGLN1-HIF1A pathway can facilitate ferroptosis. 35 We found that most genes of the prognostic System Xc-mediates the occurrence of ferroptosis by affecting glutamate uptake and glutathione synthesis. 7,36 Previous studies have shown that when glutamate is deficient, glutamate synthesis is blocked, or system Xc-is inhibited, the production of ROS and lipid peroxides decreases, which further reduces the incidence of ferroptosis. 37 Consistently, the results of functional enrichment revealed that the terms of glutamate metabolic process, glutathione peroxidase activity, regulation of iron ion transmembrane transport and immune response were mainly related to the risk model. Increasing evidence suggests that ferroptosis is associated with tumour immunity. 13,38 According to known evidence, the microenvironment of glioma is different from most solid tumours, with monocytes (macrophages and microglia) as the majority of non-neoplastic cells. 39 Glioma secretes chemokines to recruit microglia and peripherally derived monocytes, which make up as much as 30%-50% of the tumour tissue, 40 and are then differentiated into tumour-associated macrophages (TAM). TAM plays important role in promoting tumour angiogenesis, suppressing anti-tumour immune responses and assisting glioma proliferation and invasion. 41 Interestingly, we found that macrophages and monocytes had higher fractions in the highrisk group. However, there is little knowledge of the potential modulation between ferroptosis and immune regulation. Furthermore, multivariate Cox analysis showed that the risk score was an independent prognostic factor. The risk stratification of LGG is carried out according to age, surgical resection range, tumour volume, pre-operative neurological function and IDH1 status etc. 42,43 To personalize OS prediction, we combined the risk score with clinical factors to create a nomogram in LGG patients. In recent years, there have been several prediction models for glioma. 44,45 One strength of our study was that the ferroptosis-related genes of the prognostic model were closely related to the grade, TMA resistance and radiosensitization in LGG patients. This result suggests that the process of ferroptosis plays an important role in tumour differentiation and treatment sensitivity. The C-index (0.817) and calibration plots showed that our nomogram had excellent predictive power.
Therefore, physicians can apply the nomogram to improve the accuracy of identifying high-risk patients and realize accurate treatment.
In conclusion, we applied TCGA and CGGA data sets to construct and verify the prognostic signature of ferroptosis-related genes, which have the potential to predict the prognosis of LGG. In addition, we will verify the reliability of the prognostic signature in the

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.