Ferroptosis-related gene signature correlates with the tumor immune features and predicts the prognosis of glioma patients

Two public repositories, the CGGA (http://www.cgga.org.cn/↗) and the Genomic Data Commons Data Portal (GDC; https://portal.gdc.cancer.gov/↗), were searched for mRNA expression and clinical data of glioma patients to be included in the present study. In total, 934 glioma patients treated by surgical excision, living longer than 30 days from diagnosis and accessible for WHO grade information were included. Relevant transcriptome sequencing data, and clinical and molecular data were downloaded for further bioinformatics and statistical analyses. For the two RNA-seq cohorts, original data in the Fragments Per Kilobase of transcript per Million (FPKM) format were transformed to the Transcripts Per Kilobase Million (TPM) format using an algorithm described in a previous study [13,14]. The baseline clinicopathological information of glioma patients included in the present study was listed in the Table 1.Twenty-eight FRGs, including ROS proto-oncogene 1 (ROS1, receptor tyrosine kinase), acyl-CoA synthetase long chain family member 4 (ACSL4), apoptosis inducing factor mitochondria associated 2 (AIFM2), activating transcription factor 4 (ATF4), BRCA1-associated protein 1 (BAP1), BCL2 apoptosis regulator (BCL2), beclin 1 (BECN1), cyclin-dependent kinase inhibitor 1A (CDKN1A), cysteine dioxygenase type 1 (CDO1), cystathionine γ-lyase (CTH), ferritin heavy chain 1 (FTH1), ferritin light chain (FTL), glutamic-oxaloacetic transaminase 1 (GOT1), glutathione peroxidase 4 (GPX4), heat shock protein family B member 1 (HSPB1), vir-like m6A methyltransferase associated (KIAA1429, also known as VIRMA), nuclear receptor co-activator 4 (NCOA4), nuclear factor, erythroid 2 like 2 (NFE2L2), NFS1 cysteine desulfurase (NFS1), OTU deubiquitinase, ubiquitin aldehyde binding 1 (OTUB1), solute carrier family 11 member 2 (SLC11A2), solute carrier family 1 member 5 (SLC1A5), solute carrier family 40 member 1 (SLC40A1), solute carrier family 7 member 11 (SLC7A11), transferrin (TF), transferrin receptor (TFRC), tumor protein p53 (TP53), and zinc finger E-box binding homeobox 1 (ZEB1) were defined from a preceding publication [15].

Six pairs of glioma core samples and para-cancerous tissues were collected to verify the protein expression levels of HSPB1 and SLC1A5. Samples were resected from six glioma patients (two WHO grade II, one WHO grade III, and three glioblastomas (GBMs)) in the Neurosurgery Department of The Second Affiliated Hospital of Nanchang University from 2020 to 2021. These glioma core samples and para-cancerous tissues were frozen in liquid nitrogen. This research has been approved by the Medical Ethics Committee of The Second Affiliated Hospital of Nanchang University. The processes of sample collection and employ were fully in accordance with the approved guidelines. Informed consents were gained from each glioma patient.

RIPA lysis buffer with 1% PMSF was used to extract total sample protein, and protein concentrations were determined by BCA assay (Beyotime, China). A total of 15 μg total protein of each sample was added in SDS/PAGE gel and separated by electrophoresis and the separated proteins were transferred to PVDF membranes (Millipore, MA, U.S.A.). Then the PVDF membranes were blocked in 10% skim milk at 25°C for 1 h, then cropped and incubated with anti-SLC1A5 (1:1000, 20350-1-AP, Proteintech, China), anti-HSPB1 (1:1000, 18284-1-AP, Proteintech), anti-GAPDH rabbit polyclonal antibody (1:4000, 10494-1-AP, Proteintech) at 4°C overnight, respectively. The HRP-conjugated Affinipure Goat Anti-Rabbit IgG (1:4000, SA00001-2, Proteintech) were used to incubate the PVDF membranes 2 h at room temperature.

The stripes of interest proteins were developed with enhanced chemiluminescence (ECL, Thermo Fisher Scientific, 32106, U.S.A.) reagents by means of GV6000M (GelView 6000pro, China). The gray intensity of HSPB1 and SLC1A5 were calculated using the ImageJ software (National Institutes of Health, U.S.A.) and standardized to the intensity of corresponding GAPDH.

The CGGA cohort acted as the training cohort for establishing the ferroptosis-related signature (FRSig). First, we used univariate Cox regression to evaluate the prognostic roles of the 28 FRGs. Then, the expression of 22 FRGs with statistical significance (|hazard ratio (HR)| > 1.00 and P<0.05) were subjected to Lasso Cox regression, which is a classical method used to construct a high-dimension expression matrix, and a prognostic predictive FRG-signature was developed based on the CGGA cohort. Finally, a prognostic model was established based on ten FRGs; the formula was described as:

in which Coef_i indicates the coefficients, x_i is the TPM value of each FRG.

Before biological process enrichment and pathway analysis, differentially expressed genes (DEGs) between the low- and high-risk glioma subgroups were identified using the ‘limma’ package [16] under the standard of |log2(fold change) | > 1 and a P-value <0.05. A total of 1632 DEGs were identified and were included in the functional enrichment analysis using the webtool ‘Metascape’ [17], a gene annotation and analysis resource, which includes Canonical Pathways, Reactome Gene Sets, Gene Ontology (GO) Biological Processes, and the Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG pathway). Gene set enrichment analysis (GSEA) [18], a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states, was used to evaluate the significantly enriched hallmarks in the high-risk glioma patients.

We performed single-sample gene set enrichment analysis (ssGSEA) [18], an extension of GSEA, to calculate an enrichment score for each sample using a specific gene set, which quantified the level of cell infiltration or immune-related activities in the glioma TME. Immune-related gene sets used in our analysis were downloaded from the MSigDB website (https://www.gsea-msigdb.org/↗). Using these gene sets, relative enrichment scores (ESs) were calculated based on the levels of corresponding abundance of non-neoplastic cell infiltrating or immune-related activities in each glioma sample.

The prognostic nomogram model was established and validated by multivariate Cox regression method using the ‘rms’ R package. Considering the prognostic power and clinical accessibility of factors with statistical significance in the multivariate Cox regression analysis, age, WHO grade, and the FRSig risk score were incorporated in the nomogram. Calibration analysis of the prognosis predictive value of the nomogram were carried out using the ‘calibrate’ function of the ‘rms’ package. To evaluate the clinical translational potential of the nomogram, we used the decision curve analysis (DCA) to calculate the clinical benefits of different factors.

The Student’s t test and one-way ANOVA test were used to compare differences of two or three glioma subgroups, respectively, and paired t test was used to compare the protein expression difference of HSPB1 and SLC1A5 between paired glioma samples and para-cancerous tissues. In consideration of the prognostic role of the FRSig, survival curves of the prognostic analysis were determined using the Kaplan–Meier method and statistical significance was established using the two-sided log-rank test. We chose the most statistically significant cut-off point using the ‘surv-cutpoint’ function of R package ‘survminer’. We applied the univariate Cox regression model to identify the prognostic role of each FRG and calculate the HR of FRGs associated with survival. To estimate the prognostic predictive performance of the FRSig, time-dependent receiver operating characteristic curve (ROC) model analysis was carried out to calculate the area under the curve (AUC) using the ‘timeROC’ package. A two-sided statistical P-value <0.05 was considered as statistically significant in the present study.

The whole design of our study was described as Figure 1. Considering the pivotal biological functions of each ferroptosis regulator in tumorigenesis and development, we investigated the relationships between each individual ferroptosis regulator and the pathological features of gliomas, including WHO grade, isocitrate dehydrogenase (IDH) mutation status, and 1p/19q co-deletion status. The expression of each ferroptosis regulator and associated WHO grades are presented as heatmaps (Figure 2A,B). The expressions of most ferroptosis regulators are significantly associated with WHO grades. Next, we evaluated the relationship between the expression level of each individual ferroptosis regulator and IDH mutation status in lower grade gliomas (LGG, Figure 2C) and GBMs (Figure 2D), respectively. The results showed that expression of ZEB1, NCOA4, SLC11A2, TFRC, SLC1A5, HSPB1, FTH1, and FTL differed significantly between LGGs with mutated IDH and wildtype-IDH LGGs in the CGGA datasets (Figure 2C). SLC1A5, FTH1, FTL, HSPB1, NCOA4, SLC11A2, ZEB1, and TFRC genes were also differentially expressed in GBM patients with and without IDH mutation in the CGGA dataset (Figure 2D). In patients with LGG harboring IDH mutation, we observed that the expression of FTL, SLC1A5, CDKN1A, CDO1, CTH, and GOT1 was highly correlated with the status of 1p/19q co-deletion (Figure 2E).

Significant correlations between WHO grades and expression of BECN1, TP53, NFE2L2, SLC1A5, HSPB1, FTL, GOT1, NFS1, SLC11A2, SLC7A11, and TF were also confirmed by quantitative analyses using the CGGA datasets (Figure 3A). For LGGs, AIFM2, FTH1, FTL, HSPB1, SLC1A5 and TFRC were upregulated in IDH-wild LGGs and NCOA4, NFS1, OTUB1, SLC11A2 and ZEB1 were downregulated (Figure 3B). And in IDH-mutation LGGs, CDO1, CTH, FTL and SLC1A5 were overexpressed in 1p/19q non-codeletion samples and BAP1, CDKN1A, GOT1, SLC11A2, SLC7A11 and TF were downregulated in 1p/19q non-codeletion samples (Figure 3C). Additionally, AIFM2, BECN1, CDKN1A, FTH1, FTL, HSPB1, NFE2L2, ROS1, SLC1A5, SLC7A11 and TFRC were upregulated in IDH-wild GBMs, and NCOA4, OTUB1, SLC11A2, TF, TP53 and ZEB1 were downregulated in IDH-wild GBMs (Figure 3D). As the WHO grade increased, the expression of BECN1, TP53, NFE2L2, SLC1A5, HSPB1, and FTL increased, while the expression of GOT1, NFS1, SLC11A2, SLC7A11, and TF decreased.

Due to the significant associations between ferroptosis regulators and clinicopathological features of gliomas, we further investigated the prognostic values of these ferroptosis regulators using univariate Cox regression analysis in the CGGA dataset (Figure 4A). In total, 22 ferroptosis regulators were identified as associated with the overall survival (OS) of glioma patients in the CGGA cohort. The gene expression of ACSL4, BCL2, GOT1, NCOA4, NFS1, SLC11A2, SLC7A11, TF, and ZEB1 was positively associated with the OS of glioma patients, while the expression of AIFM2, ATF4, BECN1, FTH1, FTL, GPX4, HSPB1, KIAA1429, NFE2L2, SLC1A5, SLC40A1, TFRC, and TP53 were negatively correlated with OS.

Based on the 22 prognostic ferroptosis regulators identified, we developed an FRSig by performing Lasso Cox regression algorithm using the CGGA dataset (Figure 4B–D). Ten of the 22 prognostic ferroptosis regulators were selected to develop the risk signature and the expression values. The regression coefficients of the ten FRGs were used to calculate a ferroptosis-related risk score for each patient in the CGGA dataset. Glioma patients with risk scores greater than the median value were defined as high-risk glioma patients, and they showed worse OS rate and time across all glioma grades and for each single glioma grade in the CGGA dataset (Figure 4E). To further validate the prognostic value of the FRSig, risk scores were also calculated for glioma patients in the TCGA dataset using the same formula. Similarly, the FRSig could distinguish glioma patients with different clinical outcomes for all grades of glioma or for each single grade (Figure 4F).

Furthermore, we explored the mRNA expression levels of SLC1A5 and HSPB1 in the GEPIA website and the protein expression levels using clinical samples between normal brain tissues and gliomas. The results showed that SLC1A5 was up-regulated in gliomas compared with normal brain tissues in both mRNA and protein levels (Supplementary Figure S1A,B), and HSPB1 was also overexpressed in glioma samples (Supplementary Figure S1C,D).

To investigate the underlying biological functions which define the survival of glioma patients, we identified genes that were significantly differentially expressed between high- and low-risk glioma patients, and then annotated their function using GO pathway analysis for biological processes using the Metascape website. The results indicated that these DEGs were enriched in several immune-related biological processes, including the humoral immune response, production of molecular mediators of the immune response, and immunoregulatory interactions between lymphoid and non-lymphoid cells (Figure 5A–C). Furthermore, our FRSig was also associated with the development of the skeletal system, extracellular matrix organization, regulation of ion transport, and regulation of cell secretion (Figure 5A–C).

Furthermore, GSEA revealed that gene hallmarks, including allograft rejection (NES = 1.812, normalized P=0.007), IL2-STAT5 signaling (NES = 1.792, normalized P=0.002), hypoxia (NES = 1.76, normalized P=0.019), glycolysis (NES = 1.844, normalized P=0.003), DNA repair (NES = 1.754, normalized P=0.043), complement (NES = 1.749, normalized P=0.012), coagulation (NES = 1.843, normalized P<0.001), and apoptosis (NES = 1.928, normalized P=0.002), were significantly associated with high-risk glioma patients (Figure 5D).

To evaluate the prognostic role of FRSig, we investigated the associations between the FRSig and the clinical and molecular features of glioma patients. We designed a heatmap to illustrate the relationship between FRG expression, clinical and molecular features, and the risk scores in the CGGA dataset (Figure 6A). Next, violin plots were used to compare the risk score levels between glioma patients from different clinical and molecular subgroups. As shown in Figure 6B, higher risk scores were significantly associated with older age, higher WHO grades, IDH wildtype mutation status, absence of 1p/19q co-deletion, MGMT unmethylated status, and chemotherapy-treated glioma patients, while risk scores were not associated with sex or radiotherapy status.

Based on the results of enrichment analysis, we found that the FRSig was associated with the immune response process, thus we speculated that our FRSig might be associated with the cancer immune features, which would help guide the application of immunotherapy to glioma in the future. Risk scores were significantly positively associated with the immune and stromal scores, and most relevant immune cell infiltrations or immune features calculated by ssGSEA algorithm were significantly correlated with the risk score (Figure 7A). Our results indicated that glioma patients with higher risk scores correlated with higher immune scores, stromal scores, CAN percent, and tumor mutation burden (TMB) (Figure 7B). Furthermore, immune checkpoint-related genes, including C–C motif chemokine ligand 2 (CCL2), CD274 (also known as PD-L1), CD276, CD4, cytotoxic T-lymphocyte associated protein 4 (CTLA4), interleukin 1 α (IL1A), interleukin 6 (IL6), killer cell lectin like receptor B1 (CD161, also known as KLRB1), lymphocyte activating 3 (LAG3), leucine aminopeptidase 3 (LAP3), programmed cell death 1 (PD-1, also known as PDCD1), and programmed cell death 1 ligand 2 (PDCDLG2), and transforming growth factor β 1 (TGFB1), were also included in our analysis. Glioma patients with higher risk scores had higher levels of immune checkpoint mRNA expression (Figure 7C).

To assess the prognostic and predictive values of the FRSig, ROC curves were used to evaluate the predictive accuracy of OS stratifying cases by age, WHO grade, and ferroptosis-related risk score. The risk score had a higher AUC value compared with stratification by age or WHO grade for the prediction of the 1-, 3-, and 5-year OS in the CGGA dataset (Figure 8A). The ROC curves performed in the TCGA dataset also revealed that the ferroptosis-related risk score could suitably predict the 1-, 3-, or 5-year OS of gliomas with all AUC values of more than 0.800 (Figure 8B).

Moreover, we performed univariate Cox regression to evaluate the prognostic roles of clinicopathological and molecular features (including age, sex, WHO grade, IDH mutant status, 1p/19q co-deletion status, and the ferroptosis-related risk score) in both the CGGA training and TCGA validation datasets. The molecular status of IDH mutation and 1p/19q co-deletion were protective factors, while older age, WHO grade, and the ferroptosis-related risk scores were risk factors for glioma patients in the CGGA dataset (Figure 8C). We then excluded the patient sex from the multivariate Cox regression analysis given it was not significantly associated with a prognostic role in the univariate Cox regression analysis. Multivariate Cox regression showed that age, WHO grade, 1p/19q co-deletion status, and the ferroptosis-related risk score were independent prognostic factors in the CGGA dataset (Figure 8C). Similar results were found in the TCGA validation dataset, whereby age, WHO grade and FRSig were identified as independent indicators of OS of glioma patients (Figure 8D).

To further evaluate the clinical application of our FRSig, we constructed a nomogram prognostic model, based on the multivariate Cox regression results and the accessibility of clinicopathological factors in clinical diagnosis and treatment. We included the WHO grade, age, and FRSig as variables in the CGGA dataset (Figure 9A). In the CGGA cohort, calibration curves showed that the clinical nomogram model could precisely predict the 1-, 3-, and 5-year OS of patients with gliomas (C-index = 0.762) (Figure 9B), and the predictive accuracy of this nomogram was well validated in the TCGA cohort (C-index = 0.846) (Figure 9C).

Furthermore, we used DCA curves to quantify the relevant clinical benefit of the FRSig-based nomogram model. In comparison with the FRSig, age, and WHO grade, glioma patients obtained greater net benefit with a wider threshold probability for 1-, 3-, and 5-year OS prediction (Figure 9D–F) using the nomogram in the CGGA dataset. The associated DCA curves were also plotted for TCGA cohort (Figure 9G–I) and similar results were obtained, indicating that the nomogram could stably enhance the clinical benefits of glioma patients.

The original contributions presented in the study are included in the CGGA (http://www.cgga.org.cn/↗), the University of California, Santa Cruz Xena browser (UCSC Xena; https://xenabrowser.net/datapages/↗), and the GDC (https://portal.gdc.cancer.gov/↗).