Ferroptosis-Related Gene Contributes to Immunity, Stemness and Predicts Prognosis in Glioblastoma Multiforme

Ferroptosis, a recently discovered regulated programmed cell death, is associated with tumorigenesis and progression in glioblastoma. Based on widely recognized ferroptosis-related genes (FRGs), the regulation of ferroptosis patterns and corresponding characteristics of immune infiltration of 516 GBM samples with GSE13041, TCGA-GBM, and CGGA-325 were comprehensively analyzed. Here, we revealed the expression, mutations, and CNV of FRGs in GBM. We identified three distinct regulation patterns of ferroptosis and found the hub genes of immunity and stemness among DEGs in three patterns. A prognostic model was constructed based on five FRGs and verified at the mRNA and protein level. The risk score can not only predict the prognosis but also the degree of immune infiltration and ICB responsiveness by functional annotation. The overall assessment of FRGs in GBM patients will guide the direction of improved research and develop new prognostic prediction tools.

Introduction

Grade IV glioma, which is termed as glioblastoma multiforme (GBM), is the most lethal glioma (1). Despite advances in the treatment of GBM with surgery, radiation, and chemotherapy, the survival rate of patients remains 18 months (2). Previous investigations have depicted some malignant biological features that contribute to the highly recurrent and drug-resistance of GBM (3). Tremendous research studies have focused on molecular markers that contribute to GBM stemness and immunity (4, 5). Our previous study conducted a comprehensive analysis of the stemness of GBM (6). However, many therapies targeting these molecular markers become less effective in clinical practice. Therefore, the novel and effective prognostic models for the prediction of GBM prognosis and immunotherapy response need to be investigated and clarified.

Ferroptosis is a new type of programmed cell death proposed by Stockwell et al. (7). Research on the significance of ferroptosis in cancer has recently gained momentum, whereas disruption of this process under human intervention may show clinical effects (8, 9). Ferroptosis manifests cell membrane rupture and blebbing, mitochondrial and morphological changes, with the cell nucleus remain intact (10). For instance, downregulation of SLC1A5 provides melanoma cell partial immunity to ferroptosis induction (11). GOT1 inhibition promotes pancreatic cancer cell death by potentiating the activity of ferroptosis (12). Knockdown of TFRC can inhibit the cell proliferation of BRCA cell lines (13). More and more genes related to ferroptosis have been identified in glioma, such as ACSL4 that protects glioma cells and exerts antiproliferative effects by activating a ferroptosis pathway (14), ATF3 that contributes to brucine-induced glioma cell ferroptosis (15), and COPZ1 that manipulates NCOA4 to regulate the ferroptosis process in GBM (16). These ferroptosis-related genes (FRGs) are closely linked to tumorigenesis and progression. However, whether these genes are associated with the prognostic value and molecular functions of GBM patients has not been elucidated.

Despite its implication in cell death, recent studies also evaluated ferroptosis-associated diseases and their role on immunity. For example, CD8⁺ T cells suppress tumor growth by inducing ferroptosis and pyroptosis (17). In addition, ferroptosis could release various damage-associated molecular patterns (DAMPs) or lipid metabolites that are involved in the cellular immune response (18). Notably, ferroptosis was associated with tumor immune checkpoints in clear cell renal cell carcinoma (19). These researches explored the mechanisms in ferroptosis and immune microenvironment. However, these studies have not specifically focused on GBM, and the relationship between ferroptosis and immune response in GBM has not been well characterized.

Herein, we integrated data from the public Gene Expression Omnibus (GEO), The Cancer Genome Atlas (TCGA) databases, and the Chinese Glioma Genome Atlas (CGGA) to evaluate the role of FRGs signature in the prognosis in GBM patients. We further identified three distinct regulations of ferroptosis. Comparison of the DEG of three patterns unveiled five key genes involved in immunity and stemness. These genes may have potential value in the regulation of ferroptosis in GBM. Finally, a risk score based on FRGs had been constructed. Within functional annotation, we found that the risk score is not only a good predictive value for survival but also a potential factor for immune checkpoint blockade (ICB) responsiveness. Our study could help to guide the link between ferroptosis and GBM stem cells intensive research in the future and identify new ferroptosis-related targets and immune therapies.

Materials and Methods

Data Acquisition

Raw RNA-seq data (FPKM files) and clinical data on GBM were extracted from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/) and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo). After data filtration, GSE13041, TCGA-GBM, and CGGA-325 with 516 GBM tissue samples were gathered in this study for further analysis. Infiltration estimation for all TCGA tumors was collected from TIMER2.0 (http://timer.comp-genomics.org/). Copy number variant (CNV) data and somatic mutations of all genes were downloaded from the UCSC Xena browser (https://xenabrowser.net). CNV differences of all genes were calculated by the chi-square test (p < 0.05). The location of the significantly different genes on the chromosomes was shown by the RCircos R package. The protein–protein interaction network was produced by the STRING (https://www.string-db.org) database and was reconstructed via Cytoscape software. The protein expressions in human normal tissues and tumor tissues were validated via the Human Protein Altas (HPA, https://www.proteinatlas.org/).

Identification of FRGs

Ferroptosis-related genes had been categorized according to the existing literature which contains iron metabolism, oxidant metabolism, lipid metabolism, energy metabolism, and other unclassified factors (9, 20–22). According to the description of FRGs in glioma research (23), 59 genes were incorporated into follow-up studies and were provided in Supplementary Table 1. Considering the small number of normal brain tissues in the TCGA, we investigated the expression of the FRGs on the online web server GEPIA (http://gepia.cancer-pku.cn/). Differentially expressed genes (DEGs) were calculated using the R package “LIMMA” (|logFC| > 1 and p < 0.05).

Functional Enrichment Analyses

To functionally annotate DEG sets during the analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG), pathway analysis, and Gene Ontology (GO) were performed in R software version 4.0.3 using ClusterProfiler package. To calculate mRNA expression-based stemness index (mRNAsi), we used the OCLR algorithm constructed by Malta's team (24). The mRNAsi was represented using an index between zero to one to signal that the higher the mRNAsi, the greater activity of cancer stem cells. The CytoHubba plugin version 0.1 in Cytoscape version 3.8.2 was employed to identify hub genes, and enrichment analysis was performed using the ClueGO plugin version 2.5.8.

The TIMER, CIBERSORT, QUANTISEQ, Microenvironment Cell Populations-counter (MCP-counter), XCELL, and Estimating the Proportion of Immune and Cancer cells (EPIC) algorithms were used to estimate the abundance of immune cells between the high- and low-risk groups. The “ESTIMATE” R package was used to assess immune infiltration (based on the ImmuneScore, StromalScore, and ESTIMATEScore). The clustering was performed using WGCNA and the module–trait correlations with mRNAsi, EREG-mRNAsi, and ESTIMATEScore. According to the number of the genes, the minModuleSize of the mRNA was set to 50. Gene sets that could predict the responses to immune checkpoint blockade therapy were obtained from the work by Mariathasan. Single-sample gene set enrichment analysis (ssGSEA) was used to estimate immune-related functions in TCGA-GBM patients utilizing gene set variation analysis (GSVA) (25) version 1.40.1.

Tumor immune dysfunction and exclusion (TIDE) (http://tide.dfci.harvard.edu/), a well-established algorithm was employed to predict the clinical response to ICB therapy (26). TIDE is a computational framework developed to evaluate the potential of tumor immune escape from the gene expression profiles of cancer samples. The TIDE score could serve as a surrogate biomarker to predict response to ICB. The SubMap (https://www.genepattern.org/) was employed to validate the reliability of the prediction of TIDE. Mapping result is represented as a subclass association matrix filled with p-values for each subclass association (27).

Construction of a Scoring System and Calculate the Risk Score

Univariate Cox regression analysis was implemented to filtrate the prognostic FRGs. The ConsensusClusterPlus package in R was employed to investigate the detailed information in unsupervised subclasses discovery and to divide samples into appropriate parts for maximum stability (28). Thereafter, we used the R package “glmnet” to conduct least absolute shrinkage and selection operator (LASSO) Cox regression algorithm and development of a potential risk signature. The minimum value of lambda was derived from 1,000 crossvalidations (“1-se” lambda), which corresponding partial likelihood deviance value was the smallest for the risk model (29, 30). At last, coefficients with regression were confirmed by the “cvfit” function with 1,000 repeats. FRG prognostic signature involves five genes. The risk score calculating formula is as follows:

$\begin{array}{l}Riskscore={\displaystyle \sum _{i=1}^n}Coe{f}_i*x_i\end{array}$

where Coef_i means the coefficients and x_i is the expression value of each FRGs. This formula was used to calculate the risk score for each GBM patient. The predictive ability of prognostic signature for clinical traits and survival was reflected by receiver operating characteristic (ROC) and the area under the curve (AUC).

The independent clinical factor validated by univariate and multivariate Cox regression analyses was enrolled to construct a nomogram for prognosis prediction. Patients with missing data were excluded. The nomogram was performed using the “survival” and “regplot” packages of R 4.1.0 to investigate the probability of 1-, 3-, and 5-year overall survival (OS).

RNA Extraction and Real-Time PCR

For RNA extraction, three GBM tissues and one peritumoral brain edema were collected in the Second Affiliated Hospital of Harbin Medical University. This research was approved by all the patients and the ethics committee of hospital. Total RNA was isolated using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. According to the manufacturer's instructions of the Nanodrop ND-2000 spectrophotometer (Thermo Scientific, USA), 2 μg of the total RNA was transcribed into cDNA. SYBR Green PCR kit (Takara, Japan) was used for qRT-PCR. The 2–ΔΔCq method was used to calculate gene transcription level, with β-actin mRNA as control. Data represent the mean ± SD of triplicate real-time PCR. Primers (Tsingke Biotechnology Co., Ltd, Beijing, China) used are displayed in Supplementary Table 2. Clinical characteristics of patient cohort are displayed in Supplementary Table 3.

Statistical Analysis

All the data were analyzed using the R software version 4.1.0. The OS of the patients with glioma between different groups was analyzed using Kaplan–Meier curves with the log-rank test. Correlations were assessed via Spearman's coefficient. Kruskal–Wallis tests were applied for the comparison of gene expression in two or more groups. The landscape of CNV and gene location were visualized by the RCircos R package (31). A p < 0.05 was considered as statistically significant. Statistical analyses were performed using GraphPad Prism 9 for rest of the data.

Results

Landscape of FRGs in GBM

A total of 59 FRGs previously reported were included in this study. We first analyzed the expression of these genes in TCGA-GBM and normal tissues. The expression of 59 genes showed significant differences in TCGA-GBM with normal samples (Figure 1A). Among them, FANCD2, STEAP3, HMOX1, and other eight genes were upregulated in GBM (p < 0.001), whereas ACSL4, GLS2, and PEBP1 were the opposite (p < 0.001). We next examined CNV and chromosome location. Chromosome 10 carried the largest number of genes that undergo copy number variation. PGD and SLC1A5, the genes with the highest frequency of copy number loss, were located on chromosomes 1 and 19, respectively (Figures 1B,C). After that, we investigated the mutation frequencies of these genes in the TCGA-GBM dataset. As a result, there were 21 FRGs with mutation frequency >1%, and TP53 had the highest mutation frequency which was predominantly missense mutation (Figure 1D). Given the high frequency of copy number loss and mutation of TP53, we further explored the gene expression of FRGs between TP53 wildtype and mutant type. Four genes differentially expressed between subgroups were shown, and ABCC1 exhibited an increased expression in the TP53 mutant group (Figures 1E–H), which may be associated with malignant progression of TP53 mutant status (32).

FIGURE 1

Figure 1. Landscape of FRGs in TCGA-GBM. (A) The expression levels of FRGs in the TCGA dataset. Red: upregulated genes; blue: downregulated genes. (B) CNV frequency; red: CNV gain; green: CNV loss. (C) The CNV distribution of all chromosomes. (D) The mutation of FRGs in GBM; green boxes represent missense mutations, orange splicing mutations, blue frameshift mutations, and red nonsense mutations. (E–H) Four genes (ABCC1, GPX4, PEBP1, and ACSL3) with significantly different expressions between TP53 wildtype and mutant (p < 0.01).

The Relationship Between FRGs and Prognosis

Three mRNA-seq datasets that include TCGA-GBM (n = 161), CGGA325 (n = 137), and GSE12041 (n = 218) were integrated to interrogate the prognostic significance of FRGs. According to the previous studies of different metabolic pathways, genes related to ferroptosis were preliminarily divided into five categories. Gene expression, correlation, and prognostics are shown in Figure 2A. Among them, AKR1C1, AKR1C3, FDFT1 that involved in lipid metabolism, and NCOA4 that involved in iron metabolism were significantly associated with improved prognosis, which can be regarded as the protective factors. In contrast, the expression of STEAP3, HMOX1 that involved in the iron metabolism, and HSPB1 and SAT1 belong to other categories was associated with poor prognosis, which can be regarded as the risk factors (Figure 2A). Next, the Kaplan–Meier survival curve was used, and six genes that include NCOA4, STEAP3, AKR1C1, AKR1C3, FDFT1, and HSBP1 were most significantly related to OS (Figures 2B–G), which indicates that they may be vital in predicting patient prognosis.

FIGURE 2

Figure 2. FRGs associated with prognosis. (A) The degree distribution of the prognostic-related gene network. The yellow circles indicate high-risk genes, and the green circles indicate low-risk genes. The remaining five colored dots represent five types of FRGs. The blue line indicates a negative correlation, and the pink lines indicate positive correlations with correlation coefficients of p < 0.0001. Cox regression analysis (Cox p-value range 0.0001–1). The size of the dot reflects the p-value. (B–G) Six FRGs significantly associated with prognosis (NCOA4, STEAP3, AKR1C1, AKR1C3, and HSPB1 p < 0.001; FDFT1 p < 0.01).

GSVA and ssGSEA Analysis in three Clusters of FRGs

To characterize the functions of these FRGs in GBM, they were clustered for further analysis (Figure 3A). The consensus distributions for k (2 to 9) were displayed in the empirical cumulative distribution function (CDF) plots (Figures 3B,C); given the consensus matrix, k = 3 seemed to be the most suitable choice. Besides, to verify the effectiveness of unsupervised clustering, principal component analysis (PCA) can clearly show the distinction between 3 clusters which proved the accuracy of our selection for k and the effectiveness (Figure 3D). K–M analysis found significant differences in OS in the three clusters (p = 0.008), and cluster A seemed to have the poorest prognosis. Next, the expression of FRGs in the three clusters and their clinical characteristics were shown in the heatmap (Figure 3F). To gain insights into the functional implication, GSVA was performed to analyze the differentially enriched KEGG pathways in two of any three clusters. Samples in cluster A showed prominent enrichment of nod-like receptor (NLR) signaling pathway, apoptosis, amino sugar and nucleotide sugar metabolism, and cytokine–cytokine receptor interaction, etc. (Figures 3G–I). Finally, enrichment of immune cell fractions in the tumor immune microenvironment (TIME) was assessed using the ssGSEA algorithm. As a result, most types of immune cells were significantly enriched in cluster A, such as activated CD8 T cell and eosinophil (Figure 3J). These findings indicated that different regulatory patterns based on FRGs reflected the mechanisms in tumor growth, apoptosis, and immune infiltration.

FIGURE 3

Figure 3. The unsupervised clustering process for FRGs. (A) Consensus clustering matrix for k = 3. (B) Consensus clustering CDF for k = 2–9. (C) Relative change in area under the CDF curve for k = 2–9. (D) A PCA plot of unsupervised clustering when k = 3. (E) Prognostic differences between the three clusters after merging survival information (p < 0.01). (F) The heatmap of three clusters and their clinical characteristics. (G–I) Visualization for the results of KEGG for DEGs with three clusters. (J) The relative enrichment of each immune cell fraction in the TIME with the gene sets using ssGSEA (*p < 0.05, **p < 0.01, ***p < 0.001).

Identification of Hub Genes in DEGs Using WGCNA and Functional Annotation

To investigate the specific phenotype-related genes for each regulatory pattern of ferroptosis, we used the “LIMMA” package to identify the DEGs. A total of 1,622 DEGs were picked out in three clusters (Figure 4A). Functional enrichment analysis found that the main function of these DEGs enriched in neutrophil activation neutrophil-mediated immunity (BP), collagen-containing extracellular matrix, secretory granule lumen (CC), and actin-binding (MF) (Figures 4B,C). For the KEGG pathway, the most significant pathway for these DEG enrichment was cytokine–cytokine receptor interaction. Next, WGCNA was performed to structure gene coexpression networks and further identified biologically meaningful modules that corresponded to designate phenotype-related genes, which include stemness indices and the ESTIMATEScore. The most appropriate β was 7, and the relatively balanced scale independence and mean connectivity of the WGCNA were identified (Figure 4D). A total of 5 modules (merged dynamic) were identified (Figure 4E minModuleSize = 50). To analyze the correlations between merged modules and the immune and stemness phenotypes, module eigengenes (MEs), which could be regarded as representative of the gene expression patterns in a module, were determined and used to calculate the correlations with designated phenotypes. The heatmap in Figure 4F revealed the key modules (MEyellow and MEturquoise for the mRNAsi and MEturquoise and MEbrown for the ESTIMATE). The correlation graphs were plotted to select the hub genes in these modules (Figures 4G–J). Interestingly, the module turquoise had the most positive correlation with ESTIMATEScore and most negative correlation with mRNAsi. Genes with p.MMturquoise ≥ 0.8 and GS.mRNAsi ≤ −0.5 or GS.ESTIMATEScore ≥ 0.5 were screened out (Supplementary Table 4).

FIGURE 4

Figure 4. Characteristics of potential traits in different patterns of ferroptosis regulation. (A) A Venn diagram of differential genes in three patterns (p = 0.001). (B) GO analyses of differentially expressed. (C) KEGG pathway analyses of differentially expressed. (D–F) Hierarchical clustering dendrograms of identified coexpressed genes in modules. The branches of the cluster dendrogram correspond to the different gene modules. Each leaf corresponds to a gene. Each colored row represents a color-coded module, which contains a group of highly connected genes. A total of 5 modules were identified after the merger. (F) Correlations between the gene modules and target traits including mRNAsi, EREG-mRNAsi, and ESTIMATEScore. The corresponding p-value is increased in size from blue to red. (G–J) The four modules most significantly correlated with mRNAsi and ESTIMATEScore. Cor was the coefficient indices and p was Pearson's correlation.

The STRING was used to construct PPI networks. In the module turquoise, 71 genes are negatively correlated with mRNAsi, and 134 genes are positively correlated with ESTIMATEScore (GS > 0.5 and MMturquoise > 0.8). According to the confidence, the genes with significant interaction were screened, and the filtered results were imported into the Cytoscape for network visualization. The top 30 of these genes associated with mRNAsi and ESTIMATEScore with the highest combined score in STRING are exhibited (Figures 5A,B). CytoHubba and Maximal Clique Centrality (MCC) were used to explore important nodes. Top 10 MCC values were selected, and then, the intersection was taken to get the hub genes in PPI analysis (Figures 5C,D). As the result, we found that 5 genes, which include TLR4, TLR8, TNF, CD86, ITGAM, and PTPRC, were the common hub genes in the two modules. These five key genes were analyzed by KEGG and GO enrichment using the “ClueGO” and “CluePedia” plugins for Cytoscape software (Figures 5E,F p < 0.05). We found that these genes showed the enrichment of GO terms related to microglial cell activation, positive regulation of NIK/NF-kappaβ signaling and interleukin (interleukin-8 production and regulation and positive regulation of interleukin-1 beta production). Collectively, these hub genes may be the key components of the GBM immune and stemness module that contribute to immunoregulatory functions during ferroptosis.

FIGURE 5

Figure 5. Identification of hub genes from the PPI network. (A,B) PPI network obtained using STRING and showing the top 30 highest scoring genes. The bar plots of top 30 key genes of mRNAsi (A) and ESTIMATEScore (B). (C,D) Top 10 hub genes of mRNAsi (C) and ESTIMATEScore (D) explored by CytoHubba algorithm in Cytoscape software 3.8.0. and red nodes represent hub genes. (E,F) KEGG pathway and GO analysis of five overlapping genes (TLR4, TLR8, TNF, CD86, ITGAM, and PTPRC) using the “ClueGO” and “CluePedia” plugins (p < 0.05).

Cluster Analysis of DEG Levels

Previous results revealed the special performance of the DEGs in immunity and stemness of GBM. To further study, the association between these DEGs and FRGs was clustered into three categories according to their correlations in GBM. Given the consensus matrix for the analysis, k = 3 seemed to be the most suitable choice (Figures 6A–D). Similarly, these clusters also showed significant survival differences (Figure 6F p < 0.001). The differences in survival obtained by this clustering are consistent with those using FRGs (Figure 3E). The expression of DEGs in the clustering of the two methods and their clinical characteristics were shown in the heatmap (Figure 6E), and these findings were independent of clinical traits, such as age and sex. For FRGs, we found that most genes show significant differences in DEG clusters (Figure 6G). Interestingly, among the FRGs, FTH1, STEAP3, HMOX1, and 13 other genes showed an increased expression in cluster A, which had the worst prognosis.

FIGURE 6

Figure 6. Hierarchical cluster analysis of DEGs. (A–D) The consensus clustering matrix for k = 3 was determined by CDF for k = 2–9. (E) Prognostic differences between the three clusters after merging survival information (p < 0.01). (F) The heatmap of three clusters and their clinical characteristics. (G) The expression levels of FRGs in different clusters (*p < 0.05, **p < 0.01, ***p < 0.001).

Construct the FRG Prognostic Signature

Among the genes related to prognosis (Figure 2A), 16 FRGs were correlated with the OS of GBM patients. LASSO Cox analysis was performed to establish an FRG prognostic signature. In the crossvalidation process, lambda.Min was regarded as the optimal value (Figures 7A,B). A number of 5 FRGs were identified, and corresponding coefficients were calculated. A number of 516 samples were divided into train set (TCGA-GBM and GSE13041) and test set (CGGA-325), and samples were split into high- and low-risk subgroups by the median value of the risk score. Kaplan–Meier survival curves depicted that GBM patients with increased risk scores had worse clinical outcomes (Figures 7C,D, p < 0.001 in both train and test datasets). Next, we established 3- and 5-year ROC curves and found that the risk score can effectively distinguish GBM patients with different survival statuses in train set (Figure 7E, 3-year AUC = 0.706, 5-year AUC = 0.782). The risk score and survival status distributions of the train set are shown in Figure 7F. The mortality of patients increased with the increase of the risk score. The expressions of risk genes and protective genes in these 5 genes are shown in the heatmap (Figure 7G), TFRC, and STEAP3 as the risk factors increased in the high-risk score group. Conversely, NCOA4, AKR1C1, and AKR1C3 become the protective factors. The risk scores in the ferroptosis cluster and gene cluster are shown in Figures 7H,I. Univariate and multivariate Cox regression analyses show the independent prognostic value of this risk score (Supplementary Figures 1A,B). A Nomogram model was established which contained risk score, recurrent, age, and gender to assess the survival prediction in GBM patients (Supplementary Figure 1C). A Sankey diagram is used to link clustering, scoring, and survival status (Supplementary Figure 2). Most of the surviving patients belong to the low-risk group and cluster C. Finally, considering the small number of normal samples in the dataset, we validated the expression of the five key genes by qRT-PCR with the unpaired t-test in human GBM tissues. The results of qRT-PCR (Figures 8A–E) and dates in GEPIA (Figure 8F) were consistent with the expression of protein in HPA (Figures 8G–K).

FIGURE 7

Figure 7. FRG prognostic signature. (A,B) The process of building the signature. LASSO regression was performed, calculating the minimum criteria. (C,D) Kaplan–Meier curves showed that the high-risk subgroup had worse OS than the subgroup in training set (C, p < 0.001) and test set (D, p < 0.001). (E) ROC curves showed the predictive efficiency of the risk signature on the 3-year and 5-year survival rates of training set (3-year AUC = 0.706, 5-year AUC = 0.782). (F) The distributions of risk scores and the distributions of risk scores and OS status. The green and red dots indicated the alive and dead status, respectively. (G) The heatmap based on the expression of the five genes in the high- and low-risk group. (H,I) According to the formula, the different clustered samples are scored using the coefficients.

FIGURE 8

Figure 8. Validation of risk genes (A–E). (A) AKR1C1, (B) AKR1C3, (C) NCOA4, (D) STEAP3, and (E) TFRC expression in peritumoral brain edema and GBM tissues. (F) The expression of prognostic model with five genes (NCOA4, TFRC, STEAP3, AKR1C1, and AKR1C3) in GEPIA. (G–K) Validation of five genes with immunohistochemistry from the HPA database.

Enrichment Analyses of Immune-Related Functions

The enrichment scores of immune cells and corresponding immune functions and pathways with ssGSEA were quantified for the TCGA dataset. Silico approaches that include TIMER, CIBERSORT, CIBERSORT–ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC computational were employed to quantify the immune cells in high- and low-risk groups (Figure 9A). Consequently, the fraction of B cell, CD8⁺ T cell, and M2 macrophage were significantly increased in the low-risk group. NK cell and T cell regulatory (Tregs) were enriched in the high-risk group. The MHC class I scored higher in the high-risk group. The GSVA method was used to calculate the immune event scores of the high- and low-risk groups. APC coinhibition, HLA, and type I IFN response scored higher in the low-risk group (Figure 9B, p < 0.01). Next, we explored the relationship between immune checkpoint-related genes and risk score (Figure 9C). The expression of CD44, TNFRSF14, and NRP1 in the high-risk group was significantly higher than that in the low-risk group. Given this, we introduced the TIDE algorithm to assess the efficacy of FRG signatures in predicting ICB responsiveness in GBM. Submap was used to compare the prediction results (Figure 9D). As a result, different groups in train and test sets showed comparable performance in predicting the GBM response to anti-CTLA4 therapy (p < 0.05). Finally, we use GSEA to perform GO enrichment analysis on high- and low-risk GBM patients. The samples of the high-risk group were enriched in positive regulation of transcription from RNA polymerase II promoter in response to stress (GOBP), mitotic G2 M transition checkpoint (GOBP), BHLH transcription factor binding (GOMF) (Figure 9E). In the low-risk group, enriched GO terms were cell cortex region (GOCC), negative regulation of amyloid precursor protein catabolic process (GOBP), and oxidoreductase activity acting (GOMF). In conclusion, GBM patients with high- and low-risk scores had different immune-related functions (Figure 9F).

FIGURE 9

Figure 9. Enrichment analysis of immune-related functions. (A) The TIMER, CIBERSORT, CIBERSORT.ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC algorithms were applied for the immune infiltration of the high- and low-risk groups. (B) Immune score estimated by single sample GSEA (ssGSEA) in different groups. (C) Immune checkpoints in the tumor microenvironment. (D) TIDE prediction of the association of genes prognostic signature with ICB responsiveness in train and test sets. (E,F) GO analysis for high-risk group (E) and low-risk group (F) using GSEA 4.1.0.

Discussion

Malignant glioma remains a considerable threat to human health, and the prognosis of patients with GBM is dismal (33, 34). Recently, regulated cell death has gained considerable attention in cancer, especially ferroptosis (7, 35, 36). Herein, 59 FRGs in GBM were included in this study to investigate the characteristics with expression, OS, and functions. CNV exists as a genetic polymorphism in the human genome, and CNV alters tumorigenesis by deletion or amplification of a copy number of a gene (37, 38). In this study, the highest frequency of FRG PGD and SLC1A5 CNV (loss) was located on chromosomes 1 and 19, respectively. SLC1A5 expression correlated positively with immune cells, such as tumor-infiltrating B cells, CD4⁺ T in hepatocellular carcinoma, and lower-grade glioma (39). For mutations in FRGs, TP53 mutation is one of the most frequent genetic alterations in primary glioma. Previous studies have shown that TP53 polymorphism is associated with the risk of primary glioma (40). Of note, the expression of ABCC1 was higher in mutant TP53 whereas the expression of GPX4, PEBP1, and ACSL3 were higher in the wild-type TP53 group. The finding implicated TP53 mutation status was an important link in the regulation of other FRGs.

Unsupervised cluster analysis of the expression values of FRGs identified three distinct patterns in GBM. Cluster A with the worst prognosis showed high enrichment in NLR signaling pathway apoptosis, and amino sugar and nucleotide sugar metabolism. The NLR family of receptors had been recognized as the key roles of immunity and inflammation with GBM (41). Meanwhile, a variety of immune cells, such as activated CD8 T cell and eosinophil, exhibit aggregation in cluster A. Whether a direct mechanism of immune cells on ferroptosis nodes might be of physiological relevance remains elusive. Recent studies have shown that CD8 T cells may in sensitizing tumor cells toward ferroptosis (42). In addition, 1,622 DEGs were selected from three patterns. GO and KEGG pathway analysis revealed that DEG enrichment was mainly involved in the immunity biological process. Activated neutrophils are induced by the microenvironment of GBM (43). Meanwhile, some immune-associated lncRNAs in glioma were verified to be closely related to cytokine–cytokine receptor interaction (44). Functional annotation of the hub genes identified by WGCNA illuminated the potential regulatory mechanisms by which of FRGs regulate on the immune and stemness phenotypes. Recent studies suggest a possible negative regulation between stemness and immune activation (45). Glioblastoma multiforme stem cells, characterized by self-renewal and therapeutic resistance, play vital roles in GBM (6). Using the app CytoHubba in Cytoscape, we filtered 5 hub genes in both immune and stemness PPI networks. TLRs are expressed on both immune and tumor cells, which play dual roles in countering cell proliferation, migration, invasion, and glioma stem cell maintenance responses (46). TNF-α/NF-κB signaling is closely associated with glioma proliferation (47). Also, CD86 is an unfavorable prognostic biomarker in lower-grade glioma (48). It is worth noting that PTPRC, TLR8, TLR4, and TNF all exhibit functions related to IL-8 regulation. Interleukin-8 (IL-8) has been revealed as a critical regulator of central nervous system (CNS) function and development with participation in many CNS disorders including gliomas (49, 50). These are the key genes of immunity and stemness from different ferroptosis regulation patterns. This suggests that among the multiple FRGs, some genes regulate GBM stem cells and the immune microenvironment. The connection between them is also the direction of our next research.

Based on five FRGs, a prognostic model was established and validated in TCGA-GBM, GSE13041, and CGGA-325. These databases have authoritative gene expression and clinical information for GBM. The prognostic model contained five genes (NCOA4, TFRC, STEAP3, AKR1C1, and AKR1C3). The mRNA expression of genes was verified using qRT-PCR. Based on the HPA database, genes were verified at the protein levels. These FRGs affect many of the key processes involved in the tumorigenesis and progression of cancer, especially glioma. NCOA4 is a selective cargo receptor for the autophagic degradation of ferritin in glioma which is known as ferritinophagy (16, 51). The TFRC expression was higher in glioma (52), and the progression and oncogenicity of glioma were regulated by hsa-miR-144-3p/TFRC signaling (53). STEAP3 emerged as an important protein that induces mesenchymal transition and stem-like traits in glioma (54). AKR1C1 and AKR1C3 are members of the AKR superfamily which has been previously shown to be associated with oncogenic potential and proliferation capacity (55), and selective targeting of AKR1C proteins in GBM could delay the acquisition of resistance to TMZ of astroglioma cells (56). This prognostic model could predict tumor prognosis, and targeting these prognostic model genes may provide new ideas for the development of targeted treatment tools.

Our results demonstrated that a high-risk score was associated with a worse prognosis. Three-year and 5-year ROC curves indicated the 5-gene signature as a potential diagnostic factor in GBM patients. Moreover, outcome of the nomogram showed that risk score and age were associated with GBM prognosis, and it was consistent with the actual clinical situation. We investigated the correlation between high- and low-risk groups and immune cells with the CIBERSORT, CIBERSORT-ABS, QUANTISEQ, XCELL, MCPCOUNTER, and EPIC algorithms. Tregs were elevated in the high-risk group. The findings of this study are in line with those presented in previous studies, Tregs play important known roles in suppressing the immune response and maintaining immune homeostasis (57). Innovatively supporting that the abundance of nonpolarized M0 macrophages rather than M1 or M2 macrophages assembly in glioblastoma that contributed to the malignancy of tumor was proposed recently (58). Also, a recent mice study showed that increasing glioma-associated monocytes in intracranial murine GL261 leads to an increase in intratumoral and systemic myeloid-derived suppressor cells (59). In summary, regulation of ferroptosis in GBM patients may be important in controlling the inflammatory and immune responses. Research on immune checkpoints has now become a new hotspot. In this study, significant differences in the expression of immune checkpoints between high- and low-risk groups suggested that the sensitivity to immunotherapies is associated with a risk score. Moreover, our risk score may screen out potential ICB responders. This provides a new idea for in-vivo experiments of immunotherapy. However, GBM patients with low OS exhibit higher expression of markers characterizing immune response activity and T cell infiltration (60). Besides, the presence of the blood–brain barrier cannot be ignored for the nature of immunotherapy. In fact, considering that targeting these FRGs indirectly improves immunotherapy, many questions need answering.

This study still has some limitations. First, all the data used to construct and validate the prognostic model were obtained from publicly available datasets. These three GBM databases inevitably lead to the neglect of intra-tumor heterogeneity in different databases. As confirmed in the study, tumor heterogeneity has an important impact on diagnosis and treatment (61). A prospective study is needed to assess the potential application of the signature. Second, for the five key genes related to GBM stemness, which FRGs or pathways regulate them remains to be further elucidated. Third, although the survival benefits and immune-related biological processes with ferroptosis-related gene signature have been revealed through functional analysis, in-vivo and in-vitro experiments are needed to further elucidate the specific mechanism, preferably at the single-cell level in humans. Finally, we expect that this work will provide clues on immunity, stemness, and prognosis characteristics for future studies.

Conclusions

In summary, by analyzing the expression of ferroptosis-related genes in GBM, we identified three ferroptosis regulation patterns of GBM patients. Comparison of the DEG of three patterns and unveiled five key genes involved in immunity and stemness. A prognostic model based on five FRGs was built. The risk score can be a good predictor of prognosis and also predicts the degree of immune infiltration and ICB responsiveness.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found at: The Cancer Genome Atlas (TCGA) database, https://portal.gdc.cancer.gov/, TCGA-GBM, Chinese Glioma Genome Atlas (CGGA), http://www.cgga.org.cn/, CGGA-325, and National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO), https://www.ncbi.nlm.nih.gov/geo, GSE13041.

Ethics Statement

The studies involving human participants were reviewed and approved by Institutional Ethics Committee of the Second Affiliated Hospital of Harbin Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

SH, JD, and HZ conceived and designed the study and revised the manuscript. JD, FW, JJ, HJ, and XY provided analytical technical support. JD, NW, and JZ participated in the production of charts and pictures. SH supervised the study. JD designed and completed qRT-PCR experiments. All authors have read and approved the final manuscript.

Funding

This work was funded by the National Natural Science Foundation of China (no. 61575058).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors gratefully acknowledge TCGA, CGGA, and other public databases for their facilitation of data analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2022.829926/full#supplementary-material

Supplementary Figure 1. Independent prognostic value of the risk score in the train set and construction of the predictive nomogram. (A) Univariate Cox regression analysis of the signature. (B) multivariate Cox regression analysis. (C) A nomogram of the risk score for predicting 1-, 3-, and 5-year survival.

Supplementary Figure 2. Sankey diagram shows the results of our overall study. Different ferroptosis regulation patterns, three clusters of DEGs, risk scores, and survival status were linked together.

Supplementary Table 1. FRGs and classification.

Supplementary Table 2. Primer sequences for qRT-PCR.

Supplementary Table 3. Clinical characteristics of patient cohort.

Supplementary Table 4. Eligible genes in the modules.

References

1. Ostrom QT, Gittleman H, Truitt G, Boscia A, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2011–2015. Neuro-Oncology. (2018) 20:iii1–86. doi: 10.1093/neuonc/noy131

PubMed Abstract | CrossRef Full Text | Google Scholar

2. McCutcheon I, Preul M. Historical perspective on surgery and survival with glioblastoma: how far have we come? World Neurosurg. (2021) 149:148–68. doi: 10.1016/j.wneu.2021.02.047

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Olivier C, Oliver L, Lalier L, Vallette F. Drug resistance in glioblastoma: the two faces of oxidative stress. Front Mol Biosci. (2020) 7:620–77. doi: 10.3389/fmolb.2020.620677

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Chen P, Hsu W, Chang A, Tan Z, Lan Z, Zhou A, et al. Circadian regulator CLOCK recruits immune-suppressive microglia into the GBM tumor microenvironment. Cancer Discov. (2020) 10:371–81. doi: 10.1158/2159-8290.CD-19-0400

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Wang Z, Wang Y, Yang T, Xing H, Wang Y, Gao L, et al. Machine learning revealed stemness features and a novel stemness-based classification with appealing implications in discriminating the prognosis, immunotherapy and temozolomide responses of 906 glioblastoma patients. Brief Bioinformatics. (2021) 22:bbab032. doi: 10.1093/bib/bbab032

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Du J, Yan X, Mi S, Li Y, Ji H, Hou K, et al. Identification of prognostic model and biomarkers for cancer stem cell characteristics in glioblastoma by network analysis of multi-omics data and stemness indices. Front Cell Dev Biol. (2020) 8:558961. doi: 10.3389/fcell.2020.558961

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Dixon S, Lemberg K, Lamprecht M, Skouta R, Zaitsev E, Gleason C, et al. Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell. (2012) 149:1060–72. doi: 10.1016/j.cell.2012.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Xu T, Ding W, Ji X, Ao X, Liu Y, Yu W, et al. Molecular mechanisms of ferroptosis and its role in cancer therapy. J Cell Mol Med. (2019) 23:4900–12. doi: 10.1111/jcmm.14511

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hassannia B, Vandenabeele P, Vanden Berghe T. Targeting ferroptosis to iron out cancer. Cancer Cell. (2019) 35:830–49. doi: 10.1016/j.ccell.2019.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cao JY, Dixon SJ. Mechanisms of ferroptosis. Cell Mol Life Sci. (2016) 73:2195–209. doi: 10.1007/s00018-016-2194-1

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Luo M, Wu L, Zhang K, Hong W, Tian Z, Gutierrez L, et al. miR-137 regulates ferroptosis by targeting glutamine transporter SLC1A5 in melanoma. Cell Death Differ. (2018) 25:1457–72. doi: 10.1038/s41418-017-0053-8

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Kremer D, Nelson B, Lin L, Yarosz E, Halbrook C, Kerk S, et al. GOT1 inhibition promotes pancreatic cancer cell death by ferroptosis. Nat Commun. (2021) 12:4860. doi: 10.1038/s41467-021-24859-2

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Luan J, Zeng T, Zhang Q, Xia D, Cong R, Yao L, et al. A novel signature constructed by ferroptosis-associated genes (FAGs) for the prediction of prognosis in bladder urothelial carcinoma (BLCA) and associated with immune infiltration. Cancer Cell Int. (2021) 21:414. doi: 10.1186/s12935-021-02096-3

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Cheng J, Fan YQ, Liu BH, Zhou H, Wang JM, Chen QX. ACSL4 suppresses glioma cells proliferation via activating ferroptosis. Oncol Rep. (2020) 43:147–58. doi: 10.3892/or.2019.7419

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lu S, Wang XZ, He C, Wang L, Liang SP, Wang CC, et al. ATF3 contributes to brucine-triggered glioma cell ferroptosis via promotion of hydrogen peroxide and iron. Acta Pharmacol Sin. (2021) 42:1690–702. doi: 10.1038/s41401-021-00700-w

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhang Y, Kong Y, Ma Y, Ni S, Wikerholmen T, Xi K, et al. Loss of COPZ1 induces NCOA4 mediated autophagy and ferroptosis in glioblastoma cell lines. Oncogene. (2021) 40:1425–39. doi: 10.1038/s41388-020-01622-3

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Tang R, Xu J, Zhang B, Liu J, Liang C, Hua J, et al. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol. (2020) 13:110. doi: 10.1186/s13045-020-00946-7

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Proneth B, Conrad M. Ferroptosis and necroinflammation, a yet poorly explored link. Cell Death and Differentiation. (2019) 26:14–24. doi: 10.1038/s41418-018-0173-9

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Gao S, Ruan H, Liu J, Liu Y, Liu D, Tong J, et al. A novel ferroptosis-related pathway for regulating immune checkpoints in clear cell renal cell carcinoma. Front Oncol. (2021) 11:678–94. doi: 10.3389/fonc.2021.678694

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Stockwell BR, Angeli J, Bayir H. Ferroptosis: a regulated cell death nexus linking metabolism, redox biology, and disease. Cell. (2017) 171:273–85. doi: 10.1016/j.cell.2017.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Bersuker K, Hendricks JM, Li Z, Magtanong L, Ford B, Tang PH, et al. The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis. Nature. (2019) 575:688–92. doi: 10.1038/s41586-019-1705-2

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Doll S, Freitas FP, Shah R, Aldrovandi M, da Silva MC, Ingold I, et al. FSP1 is a glutathione-independent ferroptosis suppressor. Nature. (2019) 575:693–8. doi: 10.1038/s41586-019-1707-0

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wan RJ, Peng W, Xia QX, Zhou HH, Mao XY. Ferroptosis-related gene signature predicts prognosis and immunotherapy in glioma. CNS Neurosci Ther. (2021) 27:973–86. doi: 10.1111/cns.13654

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell. (2018) 173:338–54.e15. doi: 10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Hänzelmann S, Castelo R, Guinney J. GSVA gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. (2018) 24:1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Hoshida Y, Brunet JP, Tamayo P, Golub TR, Mesirov JP. Subclass mapping: identifying common subtypes in independent disease data sets. PLoS ONE. (2007) 2:e1195. doi: 10.1371/journal.pone.0001195

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang Y, Li H, Zhang W, Che Y, Bai W, Huang G. LASSO-based Cox-PH model identifies an 11-lncRNA signature for prognosis prediction in gastric cancer. Mol Med Rep. (2018) 18:5579–93. doi: 10.3892/mmr.2018.9567

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Goeman JJ. L1 penalized estimation in the Cox proportional hazards model. Biom J. (2010) 52:70–84. doi: 10.1002/bimj.200900028

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang H, Meltzer P, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinformatics. (2013) 14:244. doi: 10.1186/1471-2105-14-244

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Murnyak B, Huang LE. Association of TP53 alteration with tissue specificity and patient outcome of IDH1-mutant glioma. Cells. (2021) 10:2116. doi: 10.3390/cells10082116

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Grossman SA, Ye X, Piantadosi S, Desideri S, Nabors LB, Rosenfeld M, et al. Survival of patients with newly diagnosed glioblastoma treated with radiation and temozolomide in research studies in the United States. Clin Cancer Res. (2010) 16:2443–9. doi: 10.1158/1078-0432.CCR-09-3106

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ohgaki H, Kleihues P. Population-based studies on incidence, survival rates, and genetic alterations in astrocytic and oligodendroglial gliomas. J Neuropathol Exp Neurol. (2005) 64:479–89. doi: 10.1093/jnen/64.6.479

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Viswanathan VS, Ryan MJ, Dhruv HD, Gill S, Eichhoff OM, Seashore-Ludlow B, et al. Dependency of a therapy-resistant state of cancer cells on a lipid peroxidase pathway. Nature. (2017) 547:453–7. doi: 10.1038/nature23007

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ma S, Henson ES, Chen Y, Gibson SB. Ferroptosis is induced following siramesine and lapatinib treatment of breast cancer cells. Cell Death Dis. (2016) 7:e2307. doi: 10.1038/cddis.2016.208

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhang F, Gu W, Hurles ME, Lupski JR. Copy number variation in human health, disease, and evolution. Annu Rev Genomics Hum Genet. (2009) 10:451–81. doi: 10.1146/annurev.genom.9.081307.164217

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Falchi M, El-Sayed Moustafa JS, Takousis P, Pesce F, Bonnefond A, Andersson-Assarsson JC, et al. Low copy number of the salivary amylase gene predisposes to obesity. Nat Genet. (2014) 46:492–7. doi: 10.1038/ng.2939

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhao J, Yang Z, Tu M, Meng W, Gao H, Li MD, et al. Correlation between prognostic biomarker SLC1A5 and immune infiltrates in various types of cancers including hepatocellular carcinoma. Front Oncol. (2021) 11:608–41. doi: 10.3389/fonc.2021.608641

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang X, Sun CL, Hageman L, Smith K, Singh P, Desai S, et al. Clinical and genetic risk prediction of subsequent CNS tumors in survivors of childhood cancer: a report from the COG ALTE03N1 study. J Clin Oncol. (2017) 35:3688–96. doi: 10.1200/JCO.2017.74.7444

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Sharma N, Saxena S, Agrawal I, Singh S, Srinivasan V, Arvind S, et al. Differential expression profile of NLRs and AIM2 in glioma and implications for NLRP12 in glioblastoma. Sci Rep. (2019) 9:8480. doi: 10.1038/s41598-019-44854-4

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Wang W, Green M, Choi JE, Gijón M, Kennedy PD, Johnson JK, et al. CD8(+) T cells regulate tumour ferroptosis during cancer immunotherapy. Nature. (2019) 569:270–4. doi: 10.1038/s41586-019-1170-y

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Sippel TR, White J, Nag K, Tsvankin V, Klaassen M, Kleinschmidt-DeMasters BK, et al. Neutrophil degranulation and immunosuppression in patients with GBM: restoration of cellular immune function by targeting arginase I. Clin Cancer Res. (2011) 17:6992–7002. doi: 10.1158/1078-0432.CCR-11-1107

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Li X, Meng Y. Survival analysis of immune-related lncRNA in low-grade glioma. BMC Cancer. (2019) 19:813. doi: 10.1186/s12885-019-6032-3

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Su W, Han HH, Wang Y, Zhang B, Zhou B, Cheng Y, et al. The polycomb repressor complex 1 drives double-negative prostate cancer metastasis by coordinating stemness and immune suppression. Cancer cell. (2019) 36:139–55.e10. doi: 10.1016/j.ccell.2019.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Xun Y, Yang H, Kaminska B, You H. Toll-like receptors and toll-like receptor-targeted immunotherapy against glioma. J Hematol Oncol. (2021) 14:176. doi: 10.1186/s13045-021-01191-2

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Lei Q, Gu H, Li L, Wu T, Xie W, Li M, et al. TNIP1-mediated TNF-α/NF-κB signalling cascade sustains glioma cell proliferation. J Cell Mol Med. (2020) 24:530–8. doi: 10.1111/jcmm.14760

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Qiu H, Tian W, He Y, Li J, He C, Li Y, et al. Integrated analysis reveals prognostic value and immune correlates of CD86 expression in lower grade glioma. Front Oncol. (2021) 11:654350. doi: 10.3389/fonc.2021.654350

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Salmaggi A, Eoli M, Frigerio S, Silvani A, Gelati M, Corsini E, et al. Intracavitary VEGF, bFGF, IL-8, IL-12 levels in primary and recurrent malignant glioma. J Neurooncol. (2003) 62:297–303. doi: 10.1023/a:1023367223575

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Bonavia R, Inda MM, Vandenberg S, Cheng SY, Nagane M, Hadwiger P, et al. EGFRvIII promotes glioma angiogenesis and growth through the NF-κB, interleukin-8 pathway. Oncogene. (2012) 31:4054–66. doi: 10.1038/onc.2011.563

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Mancias JD, Wang X, Gygi SP, Harper JW, Kimmelman AC. Quantitative proteomics identifies NCOA4 as the cargo receptor mediating ferritinophagy. Nature. (2014) 509:105–9. doi: 10.1038/nature13148

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Hänninen MM, Haapasalo J, Haapasalo H, Fleming RE, Britton RS, Bacon BR, et al. Expression of iron-related genes in human brain and brain tumors. BMC Neurosci. (2009) 10:36. doi: 10.1186/1471-2202-10-36

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Ma Q, Wang X, Li J. LncRNA RP1-86C117 exacerbates the glioma progression and oncogenicity by hsa-miR-144-3p/TFRC signaling. Transl Oncol. (2021) 14:101215. doi: 10.1016/j.tranon.2021.101215

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Han M, Xu R, Wang S, Yang N, Ni S, Zhang Q, et al. Six-transmembrane epithelial antigen of prostate 3 predicts poor prognosis and promotes glioblastoma growth and invasion. Neoplasia. (2018) 20:543–54. doi: 10.1016/j.neo.2018.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Penning TM, Jin Y, Steckelbroeck S, Lanisnik Rizner T, Lewis M. Structure-function of human 3 alpha-hydroxysteroid dehydrogenases: genes and proteins. Mol Cell Endocrinol. (2004) 215:63–72. doi: 10.1016/j.mce.2003.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Le Calvé B, Rynkowski M, Le Mercier M, Bruyère C, Lonez C, Gras T, et al. Long-term in vitro treatment of human glioblastoma cells with temozolomide increases resistance in vivo through up-regulation of GLUT transporter and aldo-keto reductase enzyme AKR1C expression. Neoplasia. (2010) 12:727–39. doi: 10.1593/neo.10526

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Shao Q, Gu J, Zhou J, Wang Q, Li X, Deng Z, et al. Tissue tregs and maintenance of tissue homeostasis. Front Cell Dev Biol. (2021) 9:717903. doi: 10.3389/fcell.2021.717903

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Huang L, Wang Z, Chang Y, Wang K, Kang X, Huang R, et al. EFEMP2 indicates assembly of M0 macrophage and more malignant phenotypes of glioma. Aging. (2020) 12:8397–412. doi: 10.18632/aging.103147

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Chae M, Peterson TE, Balgeman A, Chen S, Zhang L, Renner DN, et al. Increasing glioma-associated monocytes leads to increased intratumoral and systemic myeloid-derived suppressor cells in a murine model. Neuro Oncol. (2015) 17:978–91. doi: 10.1093/neuonc/nou343

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Wang ZL, Wang Z, Li GZ, Wang QW, Jiang T. Immune cytolytic activity is associated with genetic and clinical properties of glioma. Front immunol. (2019) 10:1756. doi: 10.3389/fimmu.2019.01756

PubMed Abstract | CrossRef Full Text | Google Scholar

61. McGranahan N, Furness AJ, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science. (2016) 351:1463–9. doi: 10.1126/science.aaf1490

PubMed Abstract | CrossRef Full Text | Google Scholar