Comprehensive Analyses of Ferroptosis-Related Alterations and Their Prognostic Significance in Glioblastoma

Background: This study was designed to explore the implications of ferroptosis-related alterations in glioblastoma patients.

Method: After obtaining the data sets CGGA325, CGGA623, TCGA-GBM, and GSE83300 online, extensive analysis and mutual verification were performed using R language-based analytic technology, followed by further immunohistochemistry staining verification utilizing clinical pathological tissues.

Results: The analysis revealed a substantial difference in the expression of ferroptosis-related genes between malignant and paracancerous samples, which was compatible with immunohistochemistry staining results from clinicopathological samples. Three distinct clustering studies were run sequentially on these data. All of the findings were consistent and had a high prediction value for glioblastoma. Then, the risk score predicting model containing 23 genes (CP, EMP1, AKR1C1, FMOD, MYBPH, IFI30, SRPX2, PDLIM1, MMP19, SPOCD1, FCGBP, NAMPT, SLC11A1, S100A10, TNC, CSMD3, ATP1A2, CUX2, GALNT9, TNFAIP6, C15orf48, WSCD2, and CBLN1) on the basis of “Ferroptosis.gene.cluster” was constructed. In the subsequent correlation analysis of clinical characteristics, tumor mutation burden, HRD, neoantigen burden and chromosomal instability, mRNAsi, TIDE, and GDSC, all the results indicated that the risk score model might have a better predictive efficiency.

Conclusion: In glioblastoma, there were a large number of abnormal ferroptosis-related alterations, which were significant for the prognosis of patients. The risk score-predicting model integrating 23 genes would have a higher predictive value.

1 Introduction

Ferroptosis, a kind of controlled cell death triggered by excessive lipid peroxidation (Jiang et al., 2021), has been implicated in tumor suppression mechanisms (Chen et al., 2021) and may play a critical role in carcinogenesis and precision medicine (Shen et al., 2018; Liang et al., 2019). Additionally, it has been implicated in the control of a variety of tumor-associated signaling pathways (Wu et al., 2019; Lee et al., 2020). CD8⁺ T lymphocytes have been shown to be able to modulate tumor ferroptosis in studies in the field of tumor immunotherapy (Wang et al., 2019). Ferroptosis-related reports have been published in the fields of chemotherapy, radiation, and immunotherapy, demonstrating its distinct potential for tumor treatment (Yee et al., 2020; Liu et al., 2022; Zhao et al., 2022).

With the advancement of bioinformatics analysis and sequencing technology, it is becoming increasingly convenient for us to analyze the genomic alterations associated with certain diseases using publicly available data (Qu et al., 2016; Chen et al., 2021). Reports utilizing online data to uncover important genetic abnormalities in glioblastoma are frequent (Rutledge et al., 2013; Cimino et al., 2018) but seldom use ferroptosis. This study was created using R-based bioinformatics analytic tools and publicly available gene data from Internet sources.

2 Materials and Methods

The flow diagram of the study is provided in Supplementary Figure S1. The specific details were listed as follows.

2.1 Data Download Collection

2.1.1 TCGA-GBM Data Download

Gene expression data (https://tcga-xena-hub.s3.us-east-1.amazonaws.com/latest/TCGA.GBM.sampleMap%2FHiSeqV2.gz), genotype data (https://tcga-xena-hub.s3.us-east-1.amazonaws.com/latest/TCGA.GBM.sampleMap%2FGBM_clinicalMatrix), survival data of patient samples (https://tcga-xena-hub.s3.us-east-1.amazonaws.com/latest/survival%2FGBM_survival.txt.gz), TCGA-GBM mutant “maf” (mutation annotation format) file (https://portal.gdc.cancer.gov/files/da904cd3-79d7-4ae3-b6c0-e7127998b3e6), TCGA-GBM gene copy number data (https://gdc-hub.s3.us-east-1.amazonaws.com/latest/TCGA-GBM.gistic.tsv.gz), and the masked copy number segment file of TCGA-GBM were downloaded from the GDC database by the R package TCGAbiolinks (v 2.16.4).

2.1.2 CGGA Data Download

The data of CGGA, including mRNAseq_693, mRNAseq_325, and mRNA sequencing data (non-glioma as control), were downloaded from the following link: http://www.cgga.org.cn/download.jsp.

2.1.3 GSE83300 Data Download

The GSE83300 data were downloaded from the GEO database by the R package “GEOquery (v2.54.1).”

2.1.4 The Genomic Damage Information of TCGA Samples Was Mainly Collected From PMID: 29617664 (Knijnenburg et al., 2018)

mutLoad_nonsilent (TMB): silent mutation load per Mb. CNA_frac_altered (CNV): fraction of genome altered (fraction of bps belonging to “altered” segments), where “altered” was defined as having relative CN >0.1 or <−0.1. HRD_Score: homologous recombination deficiency score calculated from three scores (TAI + LST + HRD_LOH). HRD_TAI: number of subchromosomal regions with allelic imbalance extending to the telomere. HRD_LST: number of chromosomal breaks between adjacent regions of at least 10 Mb. HRD_LOH: the number of LOH regions of intermediate size (>15 MB but < whole chromosome in length) (Knijnenburg et al., 2018).

2.1.5 Ferroptosis-Involved Information

Ferroptosis was mainly derived from the FerrDb database (http://www.zhounan.org/ferrdb/) [(PMID: 32760210) and (PMID: 33330074)] (Liang et al., 2020; Zhuo et al., 2020). Then, the three were merged.

2.1.6 MSigDB Data Download

MSigDB data (Molecular Signatures Database) were downloaded from the following link: https://www.gsea-msigdb.org/gsea/msigdb/index.jsp.

2.1.7 Human (Gene Transfer Format) Files Download

The human “gtf” file (Homo_sapiens.GRCh38.99.gtf.gz) was downloaded from the Ensembl database, and then the symbol information was collected: (http://www.ensembl.org/info/data/ftp/index.html). The four data set samples were integrated, and the ComBat() function of the R package sva was used to remove the batch effect, and then subsequent analysis was performed.

2.2 Comprehensive Analysis of Ferroptosis-Related Genes

2.2.1 Comparison of Ferroptosis-Related Gene Diseases and Normal Expression

The “non-glioma as control” data downloaded from the CGGA database were regarded as the normal sample, and then the difference in ferroptosis-involved gene expression between cancer and normal samples in the integrated data was compared; “wilcox.test()” was used to detect significant differences, and “ggpubr (v0.4.0)” was used to achieve visualization; **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

2.2.2 Mutation and Copy Number Variation Analysis

The gene-level CNV data of TCGA-GBM were downloaded, the copy number alterations of ferroptosis-related genes were counted, and then the variation frequency was calculated, and the R package ggplot2 (v 3.3.2) was used to draw the statistical graph.

2.2.3 Circos Display and Principal Component Analysis of the Position on the Chromosome

From the human chromosome “gtf” file, the location information of the ferroptosis gene was extracted. Then, the R package RCircos (v1.2.1) was used to draw a gene circos map for position display. The function “prcomp()” of R was used to perform the principal component analysis (PCA) of cancer and normal samples. Then, R packages pca3d (v0.10.2) and rgl (v0.105.22) were adopted to draw the 3D version of the PCA.

2.3 Ferroptosis-Involved Cluster Analysis

2.3.1 Cluster Analysis of Ferroptosis-Involved Genes

Based on the expression data of ferroptosis-involved genes in cancer samples, the package ConsensusClusterPlus (v1.50.0) was used to perform the unsupervised clustering of ferroptosis genes. The clustering algorithm used was k-means. Then, combined with the overall survival (OS) data, the R packages including survival (v3.2-7) and survminer (v0.4.8) were used to perform univariate Cox analysis on all ferroptosis-involved genes, and the differences and expression correlations (Pearson coefficient) among every ferroptosis genes would be calculated. Finally, the aforementioned results would be visualized by the Cytoscape (v3.7.2).

2.3.2 Unsupervised Cluster Analysis of Samples

Based on the expression data of ferroptosis-related genes, the R package ConsensusClusterPlus (v1.50.0) was used to perform the unsupervised clustering of cancer samples. The clustering algorithm used was “pam,” and the distance used was “pearson.” Then, the R packages of survival and survminer were used to analyze the survival of the obtained subtypes, and then the Kaplan–Meier curve would be drawn.

2.4 Gene Set Variation Analysis Function Enrichment Analyses

Using the R package GSVA (v1.34.0), based on the KEGG data in MSigDB, functional enrichment analysis were performed on the samples, and then “limma” (v3.42.2) was used to retrieve the differential enrichment entries among subtypes, and the relevant threshold was set as “adj.p.value<0.05 & | logFC|> 0.3”; Finally, the R package ComplexHeatmap (v2.2.0) was used to draw the heatmap for visualization.

2.5 The Proportion and Difference of Immune-Infiltrating Cells in Different Ferroptosis Clusters Evaluated by Single Sample Gene Set Enrichment Analysis

The R package “GSVA” was used to calculate the enrichment score of 28 immune-infiltrating cells in cancer samples. After the results were obtained, the data would be normalized by the “scale()” function. According to the formula “(x-min(x))/(max(x)-min(x),” the data would be distributed from 0 to 1, and then the “wilcox.test()” was used to evaluate the significance of the difference in the proportion of immune cells among different ferroptosis cluster samples. “ggpubr” was used to achieve visualization. **** means p < 0.0001, *** means p < 0.001, **means p < 0.01, and *means p < 0.05. Finally, the Cox univariate regression analysis was performed on the proportion of immune cells, and p < 0.05 was used as the threshold to select immune cells that were significantly related to the prognosis.

2.6 Correlation Analysis Between Different Ferroptosis Clusters and Clinical Characteristics

The proportion of age, gender, chemo_therapy, and IDH1 mutations in different subtypes would be calculated, and relevant bar graphs would be drawn. Then, the “kruskal.test()” was used to test the significant difference in feature distributions among different subtypes.

2.7 Display of Ferroptosis-Related Genes in Different Ferroptosis Clusters

The expression differences of ferroptosis-related genes in different ferroptosis cluster subtypes would be counted, and then a heatmap for visualization would be drawn.

2.8 Screening of Differentially Expressed Genes in the Ferroptosis Cluster and Enrichment Analysis of DEGs

The R package “limma” was used to obtain DEGs among different subtypes. The threshold was set at “|logFC|> 1 and adj.p.val<0.05.” Afterward, the R package “clusterProfiler (v3.14.3)” was used for the functional enrichment analysis of DEGs. “p < 0.05” and “q < 0.2” were taken as the thresholds to filter the enrichment pathway, the enrichment factor would be calculated, and then the corresponding bubble chart would be drawn. The calculation formula of the enrichment factors is:

Enrichment factors = (the number of genes enriched into the pathway in the gene set)/(total number of genes in the pathway).

2.9 Ferroptosis.Gene.Cluster Obtained Based on the Cluster Analysis of DEGs

Through unsupervised clustering of samples based on DEGs, Ferroptosis.gene.cluster was obtained, and the heatmap of DEGs in different subtypes was drawn. Then, based on Ferroptosis.gene.cluster for survival analysis, the Kaplan–Meier survival curve was drawn.

2.10 Construction of “Ferroptosis Gene Signatures” Based on the DEGs of the “Ferroptosis Cluster”

First, univariate Cox regression analyses on the DEGs in different “ferroptosis clusters” were performed, and p < 0.05 was taken as the threshold to screen out genes that were significantly related to survival. Then, the R package “randomForest (v 4.614)” was used to perform random forest screening for survival significantly related genes. The related parameters used were “mtry = 2” and “ntree = 1,000.” Then, “MeanDecreaseGini> 0.72” was taken as the threshold to obtain key genes. Based on the expression of key genes, principal component analysis (PCA) on the sample was performed. The risk score is calculated by the following formula:

$\text{R}\text{S}\text{c}\text{o}\text{r}{\text{e}}_i={\displaystyle \sum \left(PC{1}_i+PC{2}_i\right)}.(1)$

Among them, PC1 and PC2 represent the scores of principal component 1 and component 2, respectively, and “i” represents the corresponding sample.

2.11 The Evaluation of the Prognostic Efficacy

After the sample risk score was obtained, the high- and low-risk score groups were divided by the median node, the survival analysis was performed in the high- and low-risk score groups, and then the Kaplan–Meier survival curve was drawn. At the same time, it was verified in the internal sub-data sets of TCGA, CGGA (CGGA325/CGGA693), and GSE83300.

After that, the time-based ROC curve was further drawn, and the AUC values of 1, 3, and 5 years were all greater than 0.6, indicating that the predictive performance of the model was better. Finally, “ggalluvial (v0.12.3)” and “ggplot2” were used to draw Sankey diagrams to express the relationship of the data characteristics.

2.12 Correlation Analyses of Risk Score, “Ferroptosis Gene,” and Pathway Function

The correlation between the risk score and the differential enrichment pathway score (Pearson correlation coefficient) and the correlation between the “ferroptosis gene” and the differential enrichment pathway score (Pearson correlation coefficient) are further calculated, and the R package “corrplot (v 0.84)” was used to complete the visualization of relevance.

2.13 Difference Analysis of Risk Score in Different Groups

2.13.1 Differences in Enrichment Scores Between High- and Low-Risk Score Groups

The differences in enrichment scores between the high- and low-risk score groups would be calculated, and the significance of the difference would be calculated by the R package “wilcox.test()” and visualized by “ggpubr.” **** means p < 0.0001, *** means p < 0.001, **means p < 0.01, and *means p < 0.05.

2.13.2 Risk Score Differences in Different Ferroptosis.gene.clusters

The risk score differences in different Ferroptosis.gene.clusters would be counted. The significance of the difference would be calculated by the R package “wilcox.test()” and visualized by “ggpubr.” **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

2.13.3 Risk Score Differences in Different “Ferroptosis Clusters”

The risk score differences in different “ferroptosis cluster” would be counted. The significance of the difference would be calculated by the R package “wilcox.test()” and visualized by “ggpubr.” **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

2.14 Risk Score Differences in Clinical Characteristics and Different Molecular Types

The distribution of risk score in age, gender, chemo_therapy, IDH1 mutation grouping, “ferroptosis cluster,” and “Ferroptosis.gene.cluster” was further checked. Then, the grouped box plot would be drawn, and the significant difference would be calculated by the R package “kruskal.test().” **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

The correlation between risk score and mutation load, homologous recombination deficiency, neoantigen load, chromosomal instability (TMB, CNV, HRD, HRD_TAI, HRD_LST, HRD_LOH, DEL, INS, and SNP), and mRNAsi would be calculated. The linear correlation graph would be drawn.

2.15 The Landscape of the High- and Low-Risk Score Groups

The R package “TCGAbiolinks (v2.16.4)” was used for the GBM’s masked copy number segment data download, and the marker file data were downloaded from the GDC Reference File (https://gdc.cancer.gov/about-data/gdc-data- processing/gdc-reference-files); Then, “GenePattern GISTIC_2.0” (https://cloud.genepattern.org/gp/pages/index.jsf) was used to analyze the alterations of CNV in the two groups online, and finally, the R package “maftools (v1.0-2)” was used for visualization.

Based on the “maf” (mutation annotation format) file of GBM mutation and risk score grouping, the mutation landscape of the two sets would be drawn by the R package “maftools.”

2.16 Immunotherapy Analysis Results in the High- and Low-Risk Score Groups (TIDE Prediction + GDSC)

2.16.1 Analysis Results of GDSC in the High- and Low-Risk Groups, Estimated by the IC₅₀ Value

The R package “pRRophetic (v0.5)” was used to predict drug treatment response in the high- and low-risk score groups, and then the box plot describing the difference would be drawn. Significant differences among groups were tested by the R package “wilcox.test().” **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

The TIDE score was used to predict the immunotherapy effect. The TIDE score could be obtained from the online website (http://tide.dfci.harvard.edu). Patients with higher TIDE scores enjoyed poorer therapeutic efficacy of immune checkpoint inhibitors and were related to the survival rate of those patients with worse anti-PD-1 and anti-CTLA-4 treatment results. After the TIDE score was obtained, the difference in different subtypes and risk groups would be calculated, and the significant difference would be tested by the R package “wilcox.test().” **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

2.16.2 Expression Differences of Immune Check Sites in the High- and Low-Risk Groups

The expression of immune checkpoints in the high- and low-risk groups would be further counted and drawn as a box-plot display; “wilcox.test()” was used to detect the significant difference. **** means p < 0.0001, *** means p < 0.001, ** means p < 0.01, and * means p < 0.05.

2.17 Immunohistochemical Validation Results of Clinical Samples

The glioblastoma samples used for this study were collected from the First Affiliated Hospital of Shandong First Medical University & Shandong Provincial Qianfoshan Hospital from June 2019 to February 2022 with informed consent provided by all participants. All tumor tissue specimens were surgically resected followed by formalin fixation and paraffin embedding (FFPE) for histological evaluation. All HE-stained and immunohistochemical (IHC)-stained slides were examined and confirmed to be glioblastoma by two experienced pathologists independently according to WHO criteria.

Slides were IHC-stained with specific primary antibodies (mouse anti-human p53 monoclonal antibody: cat. No. MAB-0674, clone MX008; mouse anti-human IDH1 R132H monoclonal antibody: cat. No. MAB-0733, clone MX031; mouse anti-human Ki-67 monoclonal antibody: cat. No. MAB-0672, clone MX006; mouse anti-human MGMT monoclonal antibody: cat. No. MAB-0361, clone MT3.1). All primary antibodies and secondary antibodies [sheep anti-mouse immunoglobulin G (IgG) polymer] were purchased from MXB Biotechnologies, Fuzhou, China. Slides were processed using an automated Roche BenchMark XT staining system according to the manufacturer’s protocol. Other genes (Ki-67, MGMT, and IDH1) were immunohistochemically stained in the same way.

3 Results

3.1 Data Collation

The flow diagram of the study is provided in Supplementary Figure S1. TCGA-GBM data were downloaded from UCSC Xena. After removing the samples without survival data, the expression matrix of 116 cancer samples was obtained. CGGA325 and CGGA693 were downloaded from the CGGA database, and then GBM data were extracted. After removing the samples without survival data, expression matrices of 137 and 237 samples were obtained, respectively. The GSE83300 data were downloaded from the GEO database; survival data and expression matrices of 50 samples were obtained. After integrating the three data sets and removing the batch effects, a total of 590 cancer samples were obtained for subsequent analyses.

Ferroptosis-related genes were downloaded from the database and two literatures (Liang et al., 2020; Zhuo et al., 2020), and finally 291 ferroptosis-related genes were obtained, of which 257 ferroptosis-related genes were displayed with expression information. The expression matrix of ferroptosis-related genes was extracted for subsequent analyses. Basic characteristics of all the data were provided in (Table 1).

TABLE 1

TABLE 1. Basic characteristics of all the data.

3.2 The Overall Display of Ferroptosis-Related Genes

3.2.1 Expression Display of Ferroptosis-Related Genes in Diseases and Normal Samples

Twenty non-glioma data were downloaded from the CGGA database and taken as the normal control, and then the expression difference of ferroptosis-related genes between the cancer and normal samples in the integrated data was compared. Among them, 239 ferroptosis-related genes were found to be significantly different between cancer and normal samples (p < 0.05), indicating that most of the expression of ferroptosis-related genes were related to GBM (Supplementary Figure S2).

3.2.2 The Mutations and CNV of Ferroptosis-Related Genes

The copy number data of ferroptosis-related genes was extracted from the gene-level copy number variation (CNV) data of TCGA-GBM (including 628 samples), and then the CNV map was drawn. Among them, CDKN2A deletion was found in more than 60% samples, and EGFR duplication was found in 40% samples (Figure 1A). The mutations of ferroptosis-related genes were further extracted from the TCGA-GBM “maf” file (containing 393 samples), and then the top 25 genes were selected and drawn into a waterfall chart, among which TP53 and EGFR ranked the top two mutation frequencies (>20%). The C > T variation in SNP was the most common (Figures 1B,C).

FIGURE 1

FIGURE 1. Mutations and copy number variations (CNVs) in ferroptosis-related genes. (A) CNV of ferroptosis-related genes. The abscissa axis represents the name of the related genes; the ordinate axis represents the CNV frequency. The type of CNV represented by red is gain; the type of CNV represented by green is loss. (B) Waterfall chart of ferroptosis-related gene mutations. The ordinate axis on the left represents the names of the top 25 genes, and the ordinate axis on the right represents the mutation frequency of the corresponding genes; different colors represent different types of gene alterations. (C) SNP of TCGA-GBM samples (C1): the abscissa axis represents the type of SNP; the ordinate axis represents the mutation percentage. (C2): the abscissa axis represents the type of variants (transitions or transversions); the ordinate axis represents the percentage of mutations. (C3): the abscissa axis represents the TCGA-GBM samples, in which different colors represent different SVP types, and the ordinate axis represents the percentage of variation in each sample.

3.2.3 Positional Circos Display on Chromosomes and Principal Component Analysis

The location information of ferroptosis-related genes was extracted from human chromosomal “gtf” files, and then a gene circos diagram was drawn for location display. Finally, the PCA results of the cancer and normal samples were drawn (Supplementary Figure S3A,B).

3.3 Ferroptosis-Related Cluster Analyses

3.3.1 Cluster Analyses of Ferroptosis-Related Genes

Based on the expression information of ferroptosis-related genes in cancer samples, the unsupervised clustering of ferroptosis-related genes was performed, and three subtypes were obtained. After that, batch Cox univariate regression analyses were performed on ferroptosis-related genes, and the Pearson coefficients among different ferroptosis-related genes were calculated at the same time. Finally, the aforementioned results were visualized by Cytoscape (Supplementary Figure S3C).

3.3.2 Unsupervised Cluster Analyses of Cancer Samples

Based on the expression information of ferroptosis-related genes, unsupervised cluster analyses of cancer samples were performed, and two cluster subtypes were obtained. After the survival analysis was performed on the two subtypes, the survival difference between the two subtypes was significant (p < 0.05) (Supplementary Figure S3D,E).

3.4 Gene Set Variation Analysis Function Enrichment Analysis

GSVA was used to perform functional enrichment analysis on the samples, and then the R package “limma” was used to retrieve differential enrichment items between the two subtypes; eight enrichment items were obtained according to the threshold, and then a heatmap was drawn. Among them, KEGG_COMPLEMENT_AND_COAGULATION_CASCADES, KEGG_ECM_ RECEPTOR_INTERACTION, KEGG_GLYCOSAMINOGLYCAN_DEGRADATION, KEGG_ GRAFT_VERSUS_HOST_DISEASE, KEGG_LEISHMANIA_INFECTION, and KEGG_OTHER_GLYCAN_DEGRADATION were found to be enriched with higher scores in cluster 1, while KEGG_PROXIMAL_TUBULE_BICARBONATE_RECLAMATION and KEGG_TERPENOID_BACKBONE_BIOSYNTHESIS were found to be enriched with higher scores in cluster 2 (Figure 2A).

FIGURE 2

FIGURE 2. (A) Enrichment score results of differential enrichment items between the two subtypes. Legend column: different colors represent different scores; cluster column: green represents cluster 1, and red represents cluster 2; type column: different colors represent the names of different databases. (B) Differences in immune-infiltrating cells in different ferroptosis clusters; the abscissa axis represents different immune-infiltrating cells; the ordinate axis represents the degree of immune cell infiltration. Cluster column: green represents cluster 1, and red represents cluster 2; ^∗∗∗∗means p < 0.0001, ^∗∗∗means p < 0.001, ^∗∗means p < 0.01, and ^∗means p < 0.05. (C) Survival risk analysis of immune-infiltrating cells. The type of immune-infiltrating cells was listed on the left side; the HR and the forest plot corresponding to immune-infiltrating cells were listed on the right side. (D) Correlation analyses between ferroptosis clusters and clinical characteristics: the distribution of different age patients in cluster 1 and cluster 2. The abscissa axis represents different cluster levels; the ordinate axis represents relative percent. Age column: red means age≤ 60; green means age >60. (E) Correlation analyses between ferroptosis clusters and clinical characteristics: the distribution of patients receiving chemotherapy in cluster 1 and cluster 2. The abscissa axis represents different cluster levels; the ordinate axis represents relative percent. Chemotherapy column: red represents no chemotherapy, and green represents chemotherapy. (F) Correlation analyses between ferroptosis clusters and clinical characteristics: the distribution of patients’ gender ratio in cluster 1 and cluster 2. The abscissa axis represents different cluster levels; the ordinate axis represents relative percent. Gender column: red represents female, and green represents male. (G) Correlation analyses between ferroptosis clusters and clinical characteristics: the distribution of IDH1 mutation status in cluster 1 and cluster 2. The abscissa axis represents different cluster levels; the ordinate axis represents relative percent. IDH1 column: red represents mutant status, and green represents wild type. (H) Correlation analyses between ferroptosis clusters and clinical characteristics: the distribution of differential enrichment pathway scores in the two subtypes. The abscissa axis represents different signal pathways; the ordinate axis represents enrichment scores. Cluster column: green represents cluster 1, and red represents cluster 2. (I) Display of ferroptosis-related genes in different ferroptosis clusters. Age column: different colors represent different age ranges; gender column: different colors represent different genders; chemotherapy column: different colors represent the status of chemotherapy; IDH1 column: different colors represent IDH1 mutation status; type column: different colors represent different data types; cluster column: different colors represent different clusters.

3.5 The Proportion and Difference of Immune-Infiltrating Cells in Different Ferroptosis Clusters Assessed by Single-Sample Gene Set Enrichment Analysis

The R package “GSVA” was used to calculate the enrichment score of 28 types of immune infiltration cells in cancer samples, and 25 types of immune infiltration cells were found to be significantly different between the two subtypes. Among them, activated CD4 T cell, activated dendritic cell, central memory CD4 T cell, effector memory CD8 T cell, gamma delta T cell, immature B cell, immature dendritic cell, macrophage, mast cell, MDSC, memory B cell, natural killer cell, natural killer T cell, neutrophil, plasmacytoid dendritic cell, regulatory T cell, T follicular helper cell, type 1 T helper cell, type 17 T helper cell, and type 2 T helper cell were found to be higher in cluster 1, while activated B cell, effector memory CD4 T cell, and monocyte were found to be higher in cluster 2 (Figure 2B).

After that, the univariate Cox regression analysis was performed on the proportion of immune cells, and the immune cells that were significantly related to the prognosis were screened with p < 0.05 as the threshold. Among them, activated dendritic cell, central memory CD4 T cell, central memory CD8 T cell, effector memory CD8 T cell, gamma delta T cell, macrophage, mast cell, MDSC, natural killer T cell, plasmacytoid dendritic cell, regulatory T cell, and T cells of both follicular helper cell and type 1 T helper cell were found to have a significant impact on the survival (Figure 2C).

3.6 The Correlation Analysis Between Different Ferroptosis Clusters and Clinical Characteristics

The proportions of age, gender, chemo_therapy, and IDH1 mutations in different subtypes would be calculated and plotted as a bar graph. Among them, only age was found to be significantly different between the two subtypes, and the remaining characteristics were not found to be significantly different. Between the two subtypes, seven of the eight differential enrichment pathway scores were found to be at a significantly different level, indicating that the two subtypes were closely related to the enrichment pathways (Figures 2D–H).

3.7 Display of Ferroptosis-Related Genes in Different “Ferroptosis Clusters”

The expression differences of ferroptosis-related genes among different “ferroptosis cluster” subtypes were counted and plotted as a heatmap (Figure 2I).

3.8 Screening for DEGs in “Ferroptosis Cluster” and Performing Enrichment Analysis for DEGs

The R package “limma” was used to screen for DEGs in different subtypes. Using |logFC|> 1 & FDR <0.05 as the thresholds, 491 DEGs were screened out, including 203 upregulated genes and 288 downregulated genes.

After that, GO and KEGG enrichment analyses of DEGs were performed. The GO enrichment analysis included three parts, namely: biological process (BP), cell component (CC), and molecular function (MF). Among them, the main pathways of BP enrichment were extracellular matrix organization and collagen fibril organization; the main pathways of CC enrichment were collagen-containing extracellular matrix and complex of collagen trimers; the main pathways of MF enrichment were extracellular matrix structural constituent, extracellular matrix structural constituent conferring tensile strength, and ion-gated channel activity. The main pathways of KEGG enrichment were ECM–receptor interaction and nicotine addiction. Relevant pathways were closely related to the ferroptosis process (Supplementary Figure S4).

3.9 “Ferroptosis.gene.cluster” Obtained by Cluster Analysis Based on DEGs

After unsupervised cluster analysis based on DEGs, two subtypes of “Ferroptosis.gene.cluster” were obtained, and the DEGs in different subtypes were displayed. Then, based on “Ferroptosis.gene.cluster” for survival analysis, the Kaplan–Meier survival curve was drawn, and the results showed that the survival curves of the two “Ferroptosis.gene.cluster” were significantly different. After that, the expression of ferroptosis-related genes in different “Ferroptosis.gene.cluster” was further analyzed, among which 196 ferroptosis-related genes were significantly different between the two “Ferroptosis.gene.clusters” (Figures 3A–C).

FIGURE 3

FIGURE 3. “Ferroptosis.gene.cluster” obtained by the clustering analysis of DEGs. (A) Heatmap of DEGs. Age column: different colors represent different age ranges; gender column: different colors represent different genders; chemotherapy column: different colors represent the status of chemotherapy; IDH1 column: different colors represent IDH1 mutation status; type column: different colors represent different data types; cluster column: different colors represent different clusters; Cluster.gene column: different colors represent different cluster.gene. (B) Kaplan–Meier survival analysis curve of Ferroptosis.gene.cluster grouping. The abscissa axis represents survival time, and the ordinate axis represents survival probability. (C) Expression differences of ferroptosis-related genes in the “Ferroptosis.gene.cluster” group. The abscissa axis represents the name of ferroptosis-related genes; the ordinate axis represents the expression level of the corresponding ferroptosis-related genes. (D) Risk score forest plot constructed by 23 key genes. The left column represents 23 key genes. The middle parts are p-value and hazard ratio. The right column is the forest plot of 23 key genes.

3.10 Construction of “Ferroptosis Gene Signatures” Based on the Risk Scores of DEGs in the “Ferroptosis Cluster”

Univariate Cox regression analysis was performed on the DEGs of different “ferroptosis clusters,” and 256 genes that were significantly related to the survival were screened out with p < 0.05 as the threshold. Then, the R package “randomForest” was used to perform random forest screening on genes that were significantly related to the survival, and 23 key genes (FMOD, MYBPH, IFI30, SRPX2, CP, PDLIM1, MMP19, TNFAIP6, SPOCD1, FCGBP, NAMPT, SLC11A1, S100A10, TNC, EMP1, C15orf48, CSMD3, GALNT9, ATP1A2, CUX2, WSCD2, AKR1C1, and CBLN1) were obtained with MeanDecreaseGini >0.72 as the threshold. Finally, based on the expression of these key genes, the PCA was performed on the samples, and the risk score of each sample was calculated (Figure 3D).

3.11 Prognostic Survival Assessment

After the sample risk score was obtained, the samples were divided into high- and low-risk score groups by the median node, and then survival analysis was performed on the two groups, and then the Kaplan–Meier curve was drawn. The significant difference of the survival could be found between the two groups. At the same time, they were verified by the CGGA (CGGA325/CGGA693), GSE83300, and TCGA internal sub-data sets. Significant differences were also found in the sub-data sets [except the TCGA data set (p = 0.056)] of the two subtype samples (Figures 4A–F).

FIGURE 4

FIGURE 4. (A–F): Kaplan–Meier survival analysis results of risk score subgroups in different data cohorts. The abscissa axis represents survival time, and the ordinate axis represents survival probability. (A) Kaplan–Meier survival analysis results of risk score subgroups in the comprehensive data cohort. (B) Kaplan–Meier survival analysis results of risk score subgroups in the CGGA data cohort. (C) Kaplan–Meier survival analysis results of risk score subgroups in the CGGA325 data cohort. (D) Kaplan–Meier survival analysis results of risk score subgroups in the CGGA693 data cohort. (E) Kaplan–Meier survival analysis results of risk score subgroups in the GSE83300 data cohort. (F) Kaplan–Meier survival analysis results of risk score subgroups in the TCGA data cohort. (G) Time-based ROC curve: the abscissa axis represents FPR (false-positive rate), and the ordinate axis represents TPR (true-positive rate). (H) Sankey diagram based on the distribution of characteristics. (I) Correlation analyses between the risk score and differential enrichment pathway score.

After that, the time-based ROC curve was drawn, and the AUC values of 1, 3, and 5 years were all greater than 0.6, indicating that the predictive model had a good predictive efficiency (Figure 4G). Finally, the Sankey diagram was drawn to indicate the relationship among data characteristics. Among them, cluster 1 of “ferroptosis clusters” had a higher proportion of deaths, and most of the survival samples had a lower risk, while the high- and low-risk samples were uniform from each sub-data set (Figure 4H).

3.12 Correlation Analysis of Risk Score, Ferroptosis-Related Genes, and Pathway Functions

The Pearson correlation coefficient between the risk score and the differential enrichment pathway scores was calculated. The results showed that most pathways were negatively correlated with the risk score, while they were positively correlated with each other. Risk score was negatively correlated with KEGG_PROXIMAL_TUBULE_BICARBONATE_RECLAMATION and KEGG_TERPENOID_BACKBONE_BIOSYNTHESIS pathways. In addition, the related pathway scores in cluster 1 were significantly lower than those of cluster 2, and the survival probability of samples in cluster 1 was lower, which was consistent with the risk score (Figure 4I).

The Pearson correlation coefficient between the ferroptosis-related genes and the differential enrichment pathway scores was further calculated and visualized by the R package “corrplot” (Figure 5A).

FIGURE 5

FIGURE 5. (A) Correlation analyses between ferroptosis-related genes and differential enrichment pathway scores. The abscissa axis represents the names of ferroptosis-related genes; the ordinate axis represents differential enrichment pathways. The legend on the right represents different Pearson correlation coefficients. (B) Difference in the enrichment scores of the subgroups with high- and low-risk scores. The abscissa axis represents pathways with different enrichment scores; the ordinate axis represents enrichment scores. Red represents the high-score group; green represents the low-score group. ^{∗∗∗∗∗}means p < 0.0001; ^∗∗∗ means p < 0.001; ^∗∗ means p < 0.01; and ^∗ means p < 0.05. (C) Risk score difference analysis between the two “Ferroptosis.gene.cluster”. The abscissa axis represents different clusters; the ordinate axis represents risk score. ^∗∗∗∗ means p < 0.0001; ^∗∗∗ means p < 0.001; ^∗∗ means p < 0.01; and ^∗ means p < 0.05. (D) Risk score difference analysis between ferroptosis clusters. The abscissa axis represents different ferroptosis clusters; the ordinate axis represents different risk scores. ^∗∗∗∗ means p < 0.0001; ^∗∗∗ means p < 0.001; ^∗∗ means p < 0.01; and ^∗ means p < 0.05. (E) Risk score difference analysis for different clinical characteristics and different molecular types. The abscissa axis represents different clinical features and molecular types; the ordinate axis represents different risk scores. ^∗∗∗∗means p < 0.0001; ^∗∗∗ means p < 0.001; ^∗∗means p < 0.01; and ^∗ means p < 0.05. (F) Difference of CNV sites in the high-risk score group. The abscissa axis represents the location of CNV on the chromosome; the ordinate axis represents G-score. (G) Difference of CNV sites in the low-risk score group. The abscissa axis represents the location of CNV on the chromosome; the ordinate axis represents G-score. (H) Mutations in the high-risk score group. The left column represents the name of the mutant genes; the right column represents the percentage of genes with mutations; different colors represent different mutation types. (I) Mutations in the low-risk score group. The left column represents the name of the mutant genes; the right column represents the percentage of genes with mutations; different colors represent different mutation types.

3.13 Differential Analysis of the Risk Score in Different Groups

3.13.1 Differential Enrichment Analysis Scores Between High- and Low-Risk Score Groups

By analyzing the differential enrichment analysis scores between the high- and low-risk score groups, it was found that eight pathways were significantly different between the high- and low-risk groups. The first six pathways in the low-risk group had significantly higher scores than those of the high-risk group. However, the scores of KEGG_PROXIMAL_TUBULE_BICARBONATE_RECLAMATION and KEGG_TERPENOID_BACKBONE_BIOSYNTHESIS in the high-risk group were significantly higher than those of the low-risk group, which was consistent with the results of the previous analysis (Figure 5B).

3.13.2 Differential Analysis of the Risk Score in Different “Ferroptosis.gene.clusters”

By analyzing the risk score differences in different “Ferroptosis.gene.clusters,” it was found that the risk score in gene.cluster 1 was significantly higher than that of gene.cluster 2 (Figure 5C), and the relevant conclusions were consistent with the previous analysis.

3.13.3 Differential Analysis of the Risk Score in Different “Ferroptosis Clusters”

Through differential analysis of the risk score in different “Ferroptosis clusters,” it was indicated that the risk score in cluster 1 was significantly higher than that of cluster 2 (Figure 5D), and the relevant conclusions were consistent with the previous analysis.

3.14 Differential Analysis of the Risk Score in Clinical Features and Different Molecular Types

The distribution of different risk scores in age, gender, chemo_therapy, IDH1 mutation grouping, “ferroptosis cluster,” and “Ferroptosis.gene.cluster” was further checked and plotted as box plots. The results showed that the distributions of different risk scores in age, IDH1 mutation group, “ferroptosis cluster,” and “Ferroptosis.gene.cluster” were significantly different (Figure 5E).

3.15 Correlation Analyses Between Risk Score and Tumor Mutation Burden, homologous recombination deficiency, Neoantigen Load, Chromosomal Instability, and mRNAsi

The correlations between risk score and TMB, HRD, neoantigen load, chromosomal instability, and mRNAsi were calculated, and then a linear correlation graph was drawn. The results indicated that risk score had a strong correlation with mRNAsi (R = −0.498), which was consistent with the previous report (Malta et al., 2018). In addition, similar significant correlation could also be found between the risk score and HRD, HRD_TAI, and HRD_LST (Supplementary Figure S5).

3.16 The Landscape of High- and Low-Risk Score Groups

The differences of CNV sites between the high- and low-risk score groups were checked, and it was found that the CNV sites of the two groups were similar. However, the G-score of 7p11.2 in the high-risk group was slightly higher than that of the low-risk group, indicating that the proportion of the relevant sites was higher in the high-risk group. In addition, 9p21.2 and 9p21.3 in the high-risk group also had higher G-scores, while 12q13.3 in the low-risk group had a higher G-score (Figures 5F,G).

Then, combined with risk score grouping, the genetic mutations of the two group samples would be displayed. The waterfall chart showed that both of the two risk groups had higher mutation rates. Among them, PTEN had the highest mutation rate in the high-risk group, while TP53 enjoyed a higher mutation rate in the low-risk group (Figures 5H,I).

3.17 Immunotherapy Analysis Results in the High- and Low-Risk Score Groups (TIDE Prediction + GDSC)

3.17.1 Analysis Results of GDSC in the High- and Low-Risk Groups, Estimated by the IC₅₀ Value

The drug treatment response of the high- and low-risk group samples was further analyzed. The results indicated that the two risk groups enjoyed significant differences in the efficacy of cisplatin, vinblastine, gemcitabine, and paclitaxel. Moreover, the samples in the low-risk group were much more sensitive to cisplatin, vinblastine, gemcitabine, and paclitaxel.

After that, the TIDE score was used to predict the immunotherapy efficacy. The results indicated that the TIDE score of the high-risk group was significantly higher than that of the low-risk group, implying that the high-risk group might have a poorer immunotherapy efficacy (Supplementary Figure S6A-E).

3.17.2 Expression Differences of Immune Check Sites in the High- and Low-Risk Groups

The expression of immune checkpoints in the high- and low-risk groups was checked. The results indicated that the expression of CD274, CTLA4, PDCD1, and LAG3 in the samples of the high-risk group was significantly higher than that of the low-risk group, suggesting that the high expression of relevant genes might be one of the reasons that affected GBM treatment efficacy (Supplementary Figure S6F–I).

3.18 Immunohistochemical Validation Results of Clinical Samples

From forty-two glioblastoma patients, including 17 females and 25 males, histopathological sections were collected. The basic characteristics of all enrolled clinical samples are summarized in Table 2. There were 15 patients with P53 gene mutation, and 27 patients were P53 wild-type. There was no significant difference in the P53 status between the two parts. IDH1 gene mutation was detected in six patients, while Ki-67 > 30% was found in 14 patients. Typical immunohistochemical staining results, including P53, Ki-67, MGMT, and IDH1, are provided in Figure 6.

TABLE 2

TABLE 2. Basic characteristics of all enrolled clinical samples.

FIGURE 6

FIGURE 6. Representative IHC analyses of p53/Ki-67/MGMT/IDH1R132H protein expression in cancer cells of glioblastoma patients. (A) Representative glioblastoma with HE staining. (B) Normal/wild-type p53 protein expression pattern with partly and weakly positive expression in tumor nuclei. Two patterns were identified as abnormal/mutant-staining pattern. (C) Abnormal overexpression of p53 protein with strong staining in nearly all tumor nuclei compared to internal control central fibroblasts. (D). Abnormal complete absence of p53 staining with sufficient staining of internal controls (fibroblasts, endothelial cells, or lymphocytes). (E) Low proportion of Ki-67 protein expression in tumor nuclei suggested that the tumor has low proliferative activity. (F) High proportion of Ki-67 protein expression in tumor nuclei suggested that the tumor has high proliferative activity. (G) Negative expression of MGMT protein in tumor nuclei might be related to MGMT methylation. (H) Strong positive expression of MGMT protein in tumor nuclei. (I) IDH1 R132H wild-type protein expression pattern with cytoplasmic negative staining of tumor cells. (J) IDH1 R132H mutation protein expression pattern with cytoplasmic positive staining of tumor cells. All mages were taken at 10 ×10 magnification on the Leica DM2000 microscope.

4.Discussion

The flow diagram of the study is provided in Supplementary Figure S1. Among the 257 ferroptosis-related genes with expression information, there were 239 ferroptosis-related genes whose expression was significantly different between cancer and normal samples (p < 0.05), indicating that the expression of most ferroptosis-related genes was related to glioblastoma (Supplementary Figure S2). The results of PCA also further verified the aforementioned conclusions (Supplementary Figure S3A,B). Some of those genes, such as CDKN2A (Kraus et al., 2001; Hu et al., 2017), CA9 (Said et al., 2007), and HSPB1 (Rajesh et al., 2020), were also found in previous reports, which further confirmed the feasibility of our study design. Further analysis found that CDKN2A, TP53, IDH1, and EGFR were still the most abnormally altered genes in glioblastoma (Figure 1) (Kraus et al., 2001; Hu et al., 2017; Hu et al., 2022; Hu et al., 2018), while the most common variation in SNP was C > T (Figure 1C) (McDonald et al., 2013; Mistry et al., 2018). This was similar to the mutation frequency of TP53 and IDH1 obtained in the immunohistochemistry of clinical samples (Table 2; Figure 6). PCA also showed that there were significant differences in ferroptosis-related genes between cancer and normal tissues (Supplementary Figure S3B). This further confirmed the important role of ferroptosis in glioblastoma.

Clustering analysis is a commonly used method for analyzing abnormal genome alterations in research (Töpfer et al., 2017; Adolfsson et al., 2019; Raymaekers and Zamar, 2020). Based on the expression information of ferroptosis-related genes in cancer samples, unsupervised clustering analysis was carried out. It was found that positive correlation was the main trend among ferroptosis-related genes (Supplementary Figure S3C), which was rarely mentioned in glioblastoma before. Through survival analysis of the two clusters obtained by the unsupervised clustering analysis of cancer samples, it was found that there was a statistically significant difference in the survival curve (Supplementary Figure S3D,E); similar differences could also be seen when GSVA, ssGSEA, and correlation analysis of clinical characteristics were completed (Figures 2A–H). Further analysis indicated that 12 types of immune cells had significant effects on survival (Figure 2C). Those would be helpful for us to use them to evaluate the clinical survival prognosis (Thomas et al., 2015; Du Four et al., 2016; Vinuesa et al., 2016; Pellegatta et al., 2018; Bayik et al., 2020; Han et al., 2020; Lupo and Matosevic, 2020; Xiong et al., 2020; Zhou et al., 2021; Campian et al., 2022).

After the differential expression analysis of ferroptosis-related genes in different ferroptosis clusters were carried out, it was found that the mutation frequency of IDH1 in cluster 2 was significantly higher than that of cluster 1 (Figures 2G–I), which further explained the reason for that cluster 2 had better survival prognosis than cluster 1 (Supplementary Figure S3E). It also further affirmed the potential prognostic significance of IDH1 mutation in glioblastoma (Nobusawa et al., 2009; Yan et al., 2009; Molenaar et al., 2014), which was also consistent with our focus on IDH1 in clinical works (Table 2; Figure 6).

Both Gene Ontology (GO) and KEGG enrichment analyses of DEGs in ferroptosis clusters suggested that ferroptosis played an important role at a certain time in the entire process of glioblastoma (Supplementary Figure S4) (Ma et al., 2019; Hynes and Naba, 2012; Bryukhovetskiy et al., 2019; Kiyokawa et al., 2021). Some of them, such as ECM–receptor interaction, nicotine addiction, and complex of collagen trimers, were rarely reported in glioblastoma (Bryukhovetskiy et al., 2019), which might provide some new ideas for further research. “Ferroptosis.gene.cluster” obtained by clustering analysis based on DEGs also showed similar prognostic effect as before (Supplementary Figure S3D,E, and Figure 3). Furthermore, the differences presented between Ferroptosis.gene.cluster 1 and Ferroptosis.gene.cluster 2 in the DEG heatmap (Figure 3A) were much more pronounced than previous clusters (Figure 2I). Statistically significant expression differences could also be found in most of the “Ferroptosis.gene.cluster” groups (Figure 3C). This means that this kind of clustering method may be a much more better approach to reveal abnormal information in glioblastoma. This was the basis for the subsequent construction of the risk score predictive model containing 23 genes. Therefore, we believed that the 23-gene risk scoring model (CP, EMP1, AKR1C1, FMOD, MYBPH, IFI30, SRPX2, PDLIM1, MMP19, SPOCD1, FCGBP, NAMPT, SLC11A1, S100A10, TNC, CSMD3, ATP1A2, CUX2, GALNT9, TNFAIP6, C15orf48, WSCD2, and CBLN1) constructed based on this might have a better prognostic prediction efficacy (Xiao et al., 2021; Yu et al., 2022), which had been verified in the subsequent analysis results (Figures 4A–G). Also, this prognostic prediction advantage was particularly evident when it was depicted in the form of “risk score” in the Sankey diagram (Figure 4H).

For those 23 genes, 3 (CP, EMP1, and AKR1C1) of them were found to be reported in ferroptosis-related studies (Yang et al., 2019; Huang et al., 2021; Huang et al., 2021), while 19 genes (CP, EMP1, AKR1C1, FMOD, MYBPH, IFI30, SRPX2, PDLIM1, MMP19, SPOCD1, FCGBP, NAMPT, SLC11A1, S100A10, TNC, CSMD3, ATP1A2, CUX2, and GALNT9) were previously reported in glioblastoma. The four tumor-related genes (TNFAIP6, C15orf48, WSCD2, and CBLN1) were neither reported in ferroptosis-related reports nor in glioblastoma-related studies (Wei et al., 2012; Su et al., 2013; Shin et al., 2020; Bushel et al., 2022), which were first found by us. While there was not necessarily a causal relationship between related things, these findings at least provided a range of options for further investigation. Therefore, the specific roles of these genes in glioblastoma still needed to be further verified by basic experiments.

We enrolled 23 genes in this risk score prediction model, which could minimize the bias in the prediction results due to the inclusion of too few predictive genes. Subsequent analysis showed that most of the pathways were negatively correlated with risk score, while these pathways were positively correlated with each other (Figure 4I). The risk score of related pathways in cluster 1 was significantly lower than that in cluster 2, and the survival probability of cluster 1 was lower, which was consistent with the results of risk score-related analysis results (Supplementary Figure S3D–3E Figure 2A, Figure 4, Figures 5A–D). This further verified the consistent trend of our overall analysis results and the feasibility of the analysis method. In the follow-up analysis results, it was found that the distribution of risk score was significantly different among age, IDH1 mutation group, “ferroptosis cluster,” and “Ferroptosis.gene.cluster” (Figure 5E). This will provide an important reference for comprehensively judging the prognosis of patients in our clinical work.

Through the correlation analysis between risk score and TMB, homologous recombination deficiency (HRD), neoantigen load and chromosomal instability, and mRNAsi (Supplementary Figure S5), we found that it had a significant correlation with HRD, HRD_TAI, and HRD_LST, especially with mRNAsi (R = −0.498), which would help us apply the mRNAsi to single-cell data to reveal patterns of intratumoral molecular heterogeneity, leading to a better understanding of glioblastoma (Malta et al., 2018; Zheng et al., 2021). By analyzing the differences of CNV loci in the high- and low-risk score groups, it was found that a G-score of 12q13.3 (mainly including OS9 gene, CDK4 gene, and SAS gene) in the low-risk group was higher (Figures 5F,G), which was considered to be the characteristic locus (Reifenberger et al., 1994; Rollbrocker et al., 1996). Regardless of grouping, missense mutations, including PTEN, TP53, EGFR, and TTN, consistent with former reports and clinical verification (Table 2 and Figure 6), were still the predominant genomic alteration type in glioblastoma (Figures 5H,I) (Smith et al., 2001; Ohgaki et al., 2004; Lee et al., 2006; Binder et al., 2018). Among them, PTEN had the highest alteration rate in the high-risk score group, while TP53 had a higher alteration rate in the low-risk score group (Figures 5H,I).

TIDE and GDSC assays were often used to assess anti-tumor therapy responses (Wang et al., 2021; Cascio et al., 2021; Cancer Cell Line Encyclopedia Consortium, 2015). In this study, the results of these two analyses showed significant differences in different risk score subgroups (Supplementary Figure S6A–E), further supporting the feasibility of the risk score model to be used for the clinical prognostic assessment. Subsequent immune checkpoint analyses revealed that the expressions of CD274, CTLA4, PDCD1, and LAG3 in the high-risk group were significantly higher than those of the low-risk group, suggesting that the high expression of related genes might be one of the reasons that affected GBM treatment and prognosis (Supplementary Figure 6F–I), which was also consistent with previous reports on the correlation between those four genes and the glioblastoma prognosis (Andrews et al., 2020; Du et al., 2020; Yang et al., 2021; Bi et al., 2022; Preddy et al., 2022).

There were three groups of consistent clustering analyses in this study. The first was the gene clustering analysis for Ferroptosis gene; the second was the clustering analysis based on the expression of Ferroptosis gene, named “Ferroptosis.cluster”; the last was the expression of DEGs based on Ferroptosis.cluster named “Ferroptosis.gene.cluster.” The whole process of cluster analysis was also an optimization process of all the data. Therefore, the consistency of the obtained results had a better feasibility and reliability. In addition, the data size of this study was relatively sufficient, and the analysis results of the verification data were consistent from each other and had a statistical significance, which might provide a certain reference for later basic or clinical research in this field.

5 Conclusion

In glioblastoma, there were a large number of abnormal ferroptosis-related alterations, which were significant for the prognosis of patients. The risk score-predicting model including 23 ferroptosis-related genes might have a better predictive significance.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

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

Author Contributions

The study was designed by YT and YS; the collection of clinical information was completed by TW; HL was responsible for data collection and immunohistochemical staining of clinical samples; all the data cleaning and analyses were put into practice by YT, WL, and XY; the manuscript was drafted by YT; and YT and YS reviewed the manuscript for scientific soundness. All authors reviewed the final draft and approved its submission.

Funding

The sources of funding for this study are listed as follows: the Academic Promotion Program of Shandong First Medical University (2019QL025; YS), the Natural Science Foundation of Shandong Province (ZR2019MH042; YS), the Jinan Science and Technology Program (201805064; YS), and the Postdoctoral Innovation Project of Jinan (YT).

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.

Supplementary Material

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

Abbreviations

TMB, tumor mutational burden; PCA, principal component analysis; GSVA, gene set variation analysis; ssGSEA, single-sample gene set enrichment analysis; FPR, false-positive rate; TPR, true-positive rate; HRD, homologous recombination deficiency; TAI, telomeric–allelic imbalance; LST, large-scale state transitions; LOH, loss of heterozygosity; CNV, copy number variation; TIDE, tumor immune dysfunction and exclusion; GDSC, Genomics of Drug Sensitivity in Cancer; MSigDB, Molecular Signatures Database.

References

Adolfsson, A., Ackerman, M., and Brownstein, N. C. (2019). To Cluster, or Not to Cluster: An Analysis of Clusterability Methods. Pattern Recognit. 88, 13–26. doi:10.1016/j.patcog.2018.10.026

CrossRef Full Text | Google Scholar

Andrews, L. P., Somasundaram, A., Moskovitz, J. M., Szymczak-Workman, A. L., Liu, C., Cillo, A. R., et al. (2020). Resistance to PD1 Blockade in the Absence of Metalloprotease-Mediated LAG3 Shedding. Sci. Immunol. 5 (49), eabc2728. doi:10.1126/sciimmunol.abc2728

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayik, D., Zhou, Y., Park, C., Hong, C., Vail, D., Silver, D. J., et al. (2020). Myeloid-Derived Suppressor Cell Subsets Drive Glioblastoma Growth in a Sex-specific Manner. Cancer Discov. 10 (8), 1210–1225. doi:10.1158/2159-8290.CD-19-1355

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, Y., Wu, Z.-H., and Cao, F. (2022). Prognostic Value and Immune Relevancy of a Combined Autophagy-, Apoptosis- and Necrosis-Related Gene Signature in Glioblastoma. BMC Cancer 22 (1), 233. doi:10.1186/s12885-022-09328-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Binder, Z. A., Thorne, A. H., Bakas, S., Wileyto, E. P., Bilello, M., Akbari, H., et al. (2018). Epidermal Growth Factor Receptor Extracellular Domain Mutations in Glioblastoma Present Opportunities for Clinical Imaging and Therapeutic Development. Cancer Cell 34 (1), 163–177. doi:10.1016/j.ccell.2018.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Bryukhovetskiy, I. S., Khotimchenko, Y., and Shevchenko, V. (2019). Molecular Determinants of Interaction between Glioblastoma CD133+ Cancer Stem Cells and Extracellular Matrix. Ann. Oncol. 30, vii4. doi:10.1093/annonc/mdz413.016

CrossRef Full Text | Google Scholar

Bushel, P. R., Ward, J., Burkholder, A., Li, J., and Anchang, B. (2022). Mitochondrial-nuclear Epistasis Underlying Phenotypic Variation in Breast Cancer Pathology. Sci. Rep. 12 (1), 1393. doi:10.1038/s41598-022-05148-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Campian, J. L., Ghosh, S., Kapoor, V., Yan, R., Thotala, S., Jash, A., et al. (2022). Long-Acting Recombinant Human Interleukin-7, NT-I7, Increases Cytotoxic CD8 T Cells and Enhances Survival in Mouse Glioma Models. Clin. Cancer Res. 28, 1229–1239. doi:10.1158/1078-0432.CCR-21-0947

PubMed Abstract | CrossRef Full Text | Google Scholar

Cancer Cell Line Encyclopedia Consortium (2015). Pharmacogenomic Agreement between Two Cancer Cell Line Data Sets. Nature 528 (7580), 84–87. doi:10.1038/nature15736

PubMed Abstract | CrossRef Full Text | Google Scholar

Cascio, S., Chandler, C., Zhang, L., Sinno, S., Gao, B., Onkar, S., et al. (2021). Cancer-associated MSC Drive Tumor Immune Exclusion and Resistance to Immunotherapy, Which Can Be Overcome by Hedgehog Inhibition. Sci. Adv. 7 (46), eabi5790. doi:10.1126/sciadv.abi5790

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Wang, C., Sun, H., Wang, J., Liang, Y., Wang, Y., et al. (2021). The Bioinformatics Toolbox for circRNA Discovery and Analysis. Brief. Bioinform 22 (2), 1706–1728. doi:10.1093/bib/bbaa001

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Kang, R., Kroemer, G., and Tang, D. (2021). Broadening Horizons: the Role of Ferroptosis in Cancer. Nat. Rev. Clin. Oncol. 18 (5), 280–296. doi:10.1038/s41571-020-00462-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cimino, P. J., McFerrin, L., Wirsching, H.-G., Arora, S., Bolouri, H., Rabadan, R., et al. (2018). Copy Number Profiling across Glioblastoma Populations Has Implications for Clinical Trial Design. Neuro Oncol. 20 (10), 1368–1373. doi:10.1093/neuonc/noy108

PubMed Abstract | CrossRef Full Text | Google Scholar

Du Four, S., Maenhout, S. K., Benteyn, D., De Keersmaecker, B., Duerinck, J., Thielemans, K., et al. (2016). Disease Progression in Recurrent Glioblastoma Patients Treated with the VEGFR Inhibitor Axitinib Is Associated with Increased Regulatory T Cell Numbers and T Cell Exhaustion. Cancer Immunol. Immunother. 65 (6), 727–740. doi:10.1007/s00262-016-1836-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, L., Lee, J.-H., Jiang, H., Wang, C., Wang, S., Zheng, Z., et al. (2020). β-Catenin Induces Transcriptional Expression of PD-L1 to Promote Glioblastoma Immune Evasion. J. Exp. Med. 217 (11), e20191115. doi:10.1084/jem.20191115

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, J., Khatwani, N., Searles, T. G., Turk, M. J., and Angeles, C. V. (2020). Memory CD8+ T Cell Responses to Cancer. Seminars Immunol. 49, 101435. doi:10.1016/j.smim.2020.101435

CrossRef Full Text | Google Scholar

Hu, C., Leche, C. A., Kiyatkin, A., Yu, Z., Stayrook, S. E., Ferguson, K. M., et al. (2022). Glioblastoma Mutations Alter EGFR Dimer Structure to Prevent Ligand Bias. Nature 602 (7897), 518–522. doi:10.1038/s41586-021-04393-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, H., Mu, Q., Bao, Z., Chen, Y., Liu, Y., Chen, J., et al. (2018). Mutational Landscape of Secondary Glioblastoma Guides MET-Targeted Trial in Brain Tumor. Cell 175 (6), 1665–1678. doi:10.1016/j.cell.2018.09.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, L. S., Ning, S., Eschbacher, J. M., Baxter, L. C., Gaw, N., Ranjbar, S., et al. (2017). Radiogenomics to Characterize Regional Genetic Heterogeneity in Glioblastoma. Neuro Oncol. 19 (1), 128–137. doi:10.1093/neuonc/now135

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, F., Zheng, Y., Li, X., Luo, H., and Luo, L. (2021). Ferroptosis-related Gene AKR1C1 Predicts the Prognosis of Non-small Cell Lung Cancer. Cancer Cell Int. 21 (1), 567. doi:10.1186/s12935-021-02267-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hynes, R. O., and Naba, A. (2012). Overview of the Matrisome-Aan Inventory of Extracellular Matrix Constituents and Functions. Cold Spring Harb. Perspect. Biol. 4 (1), a004903. doi:10.1101/cshperspect.a004903

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, X., Stockwell, B. R., and Conrad, M. (2021). Ferroptosis: Mechanisms, Biology and Role in Disease. Nat. Rev. Mol. Cell Biol. 22 (4), 266–282. doi:10.1038/s41580-020-00324-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiyokawa, J., Kawamura, Y., Ghouse, S. M., Acar, S., Barçın, E., Martínez-Quintanilla, J., et al. (2021). Modification of Extracellular Matrix Enhances Oncolytic Adenovirus Immunotherapy in Glioblastoma. Clin. Cancer Res. 27 (3), 889–902. doi:10.1158/1078-0432.CCR-20-2400

PubMed Abstract | CrossRef Full Text | Google Scholar

Knijnenburg, T. A., Wang, L., Zimmermann, M. T., Chambwe, N., Gao, G. F., Cherniack, A. D., et al. (2018). Genomic and Molecular Landscape of DNA Damage Repair Deficiency across the Cancer Genome Atlas. Cell Rep. 23 (1), 239. doi:10.1016/j.celrep.2018.03.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Kraus, J. A., Lamszus, K., Glesmann, N., Beck, M., Wolter, M., Sabel, M., et al. (2001). Molecular Genetic Alterations in Glioblastomas with Oligodendroglial Component. Acta Neuropathol. 101 (4), 311–320. doi:10.1007/s004010000258

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H., Zandkarimi, F., Zhang, Y., Meena, J. K., Kim, J., Zhuang, L., et al. (2020). Energy-stress-mediated AMPK Activation Inhibits Ferroptosis. Nat. Cell Biol. 22 (2), 225–234. doi:10.1038/s41556-020-0461-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. C., Vivanco, I., Beroukhim, R., Huang, J. H. Y., Feng, W. L., DeBiasi, R. M., et al. (2006). Epidermal Growth Factor Receptor Activation in Glioblastoma through Novel Missense Mutations in the Extracellular Domain. PLoS Med. 3 (12), e485. doi:10.1371/journal.pmed.0030485

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, C., Zhang, X., Yang, M., and Dong, X. (2019). Recent Progress in Ferroptosis Inducers for Cancer Therapy. Adv. Mat. 31 (51), 1904197. doi:10.1002/adma.201904197

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, J.-y., Wang, D.-s., Lin, H.-c., Chen, X.-x., Yang, H., Zheng, Y., et al. (2020). A Novel Ferroptosis-Related Gene Signature for Overall Survival Prediction in Patients with Hepatocellular Carcinoma. Int. J. Biol. Sci. 16 (13), 2430–2441. doi:10.7150/ijbs.45050

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Zhu, C., Chen, X., Guan, G., Zou, C., Shen, S., et al. (2022). Ferroptosis, as the Most Enriched Programmed Cell Death Process in Glioma, Induces Immunosuppression and Immunotherapy Resistance. Neuro Oncol. noac033. doi:10.1093/neuonc/noac033

PubMed Abstract | CrossRef Full Text | Google Scholar

Lupo, K. B., and Matosevic, S. (2020). CD155 Immunoregulation as a Target for Natural Killer Cell Immunotherapy in Glioblastoma. J. Hematol. Oncol. 13 (1), 76. doi:10.1186/s13045-020-00913-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, D., Liu, S., Lal, B., Wei, S., Wang, S., Zhan, D., et al. (2019). Extracellular Matrix Protein Tenascin C Increases Phagocytosis Mediated by CD47 Loss of Function in Glioblastoma. Cancer Res. 79 (10), 2697–2708. doi:10.1158/0008-5472.CAN-18-3125

PubMed Abstract | CrossRef Full Text | Google Scholar

Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., et al. (2018). Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell 173 (2), 338–354. doi:10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

McDonald, K. L., Rapkins, R. W., Olivier, J., Zhao, L., Nozue, K., Lu, D., et al. (2013). The T Genotype of the MGMT C>T (Rs16906252) Enhancer Single-Nucleotide Polymorphism (SNP) Is Associated with Promoter Methylation and Longer Survival in Glioblastoma Patients. Eur. J. Cancer 49 (2), 360–368. doi:10.1016/j.ejca.2012.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Mistry, A. M., Vnencak-Jones, C. L., and Mobley, B. C. (2018). Clinical Prognostic Value of the Isocitrate Dehydrogenase 1 Single-Nucleotide Polymorphism Rs11554137 in Glioblastoma. J. Neurooncol 138 (2), 307–313. doi:10.1007/s11060-018-2796-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Molenaar, R. J., Verbaan, D., Lamba, S., Zanon, C., Jeuken, J. W. M., Boots-Sprenger, S. H. E., et al. (2014). The Combination of IDH1 Mutations and MGMT Methylation Status Predicts Survival in Glioblastoma Better Than Either IDH1 or MGMT Alone. Neuro Oncol. 16 (9), 1263–1273. doi:10.1093/neuonc/nou005

PubMed Abstract | CrossRef Full Text | Google Scholar

Nobusawa, S., Watanabe, T., Kleihues, P., and Ohgaki, H. (2009). IDH1 Mutations as Molecular Signature and Predictive Factor of Secondary Glioblastomas. Clin. Cancer Res. 15 (19), 6002–6007. doi:10.1158/1078-0432.CCR-09-0715

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohgaki, H., Dessen, P., Jourde, B., Horstmann, S., Nishikawa, T., Di Patre, P.-L., et al. (2004). Genetic Pathways to Glioblastoma: A Population-Based Study. Cancer Res. 64 (19), 6892–6899. doi:10.1158/0008-5472.CAN-04-1337

PubMed Abstract | CrossRef Full Text | Google Scholar

Pellegatta, S., Eoli, M., Cuccarini, V., Anghileri, E., Pollo, B., Pessina, S., et al. (2018). Survival Gain in Glioblastoma Patients Treated with Dendritic Cell Immunotherapy Is Associated with Increased NK but Not CD8+ T Cell Activation in the Presence of Adjuvant Temozolomide. Oncoimmunology 7 (4), e1412901. doi:10.1080/2162402X.2017.1412901

PubMed Abstract | CrossRef Full Text | Google Scholar

Preddy, I., Nandoliya, K., Miska, J., and Ahmed, A. U. (2022). Checkpoint: Inspecting the Barriers in Glioblastoma Immunotherapies. Seminars Cancer Biol. doi:10.1016/j.semcancer.2022.02.012

CrossRef Full Text | Google Scholar

Qu, K., Garamszegi, S., Wu, F., Thorvaldsdottir, H., Liefeld, T., Ocana, M., et al. (2016). Integrative Genomic Analysis by Interoperation of Bioinformatics Tools in GenomeSpace. Nat. Methods 13 (3), 245–247. doi:10.1038/nmeth.3732

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajesh, Y., Biswas, A., Kumar, U., Banerjee, I., Das, S., Maji, S., et al. (2020). Lumefantrine, an Antimalarial Drug, Reverses Radiation and Temozolomide Resistance in Glioblastoma. Proc. Natl. Acad. Sci. U.S.A. 117 (22), 12324–12331. doi:10.1073/pnas.1921531117

PubMed Abstract | CrossRef Full Text | Google Scholar

Raymaekers, J., and Zamar, R. H. (2020). Pooled Variable Scaling for Cluster Analysis. Bioinformatics 36 (12), 3849–3855. doi:10.1093/bioinformatics/btaa243

PubMed Abstract | CrossRef Full Text | Google Scholar

Reifenberger, G., Reifenberger, J., Ichimura, K., Meltzer, P. S., and Collins, V. P. (1994). Amplification of Multiple Genes from Chromosomal Region 12q13-14 in Human Malignant Gliomas: Preliminary Mapping of the Amplicons Shows Preferential Involvement of CDK4, SAS, and MDM2. Cancer Res. 54 (16), 4299–4303.

PubMed Abstract | Google Scholar

Rollbrocker, B., Waha, A., Louis, D. N., Wiestler, O. D., and von Deimling, A. (1996). Amplification of the Cyclin-dependent Kinase 4 ( CDK4 ) Gene Is Associated with High Cdk4 Protein Levels in Glioblastoma Multiforme. Acta Neuropathol. 92 (1), 70–74. doi:10.1007/s004010050491

PubMed Abstract | CrossRef Full Text | Google Scholar

Rutledge, W. C., Kong, J., Gao, J., Gutman, D. A., Cooper, L. A. D., Appin, C., et al. (2013). Tumor-infiltrating Lymphocytes in Glioblastoma Are Associated with Specific Genomic Alterations and Related to Transcriptional Class. Clin. Cancer Res. 19 (18), 4951–4960. doi:10.1158/1078-0432.CCR-13-0551

PubMed Abstract | CrossRef Full Text | Google Scholar

Said, H. M., Hagemann, C., Staab, A., Stojic, J., Kühnel, S., Vince, G. H., et al. (2007). Expression Patterns of the Hypoxia-Related Genes Osteopontin, CA9, Erythropoietin, VEGF and HIF-1α in Human Glioma In Vitro and In Vivo. Radiotherapy Oncol. 83 (3), 398–405. doi:10.1016/j.radonc.2007.05.003

CrossRef Full Text | Google Scholar

Shen, Z., Song, J., Yung, B. C., Zhou, Z., Wu, A., and Chen, X. (2018). Emerging Strategies of Cancer Therapy Based on Ferroptosis. Adv. Mat. 30 (12), 1704007. doi:10.1002/adma.201704007

PubMed Abstract | CrossRef Full Text | Google Scholar

Shin, S.-B., Jang, H.-R., Xu, R., Won, J.-Y., and Yim, H. (2020). Active PLK1-Driven Metastasis Is Amplified by TGF-βSignaling that Forms a Positive Feedback Loop in Non-small Cell Lung Cancer. Oncogene 39 (4), 767–785. doi:10.1038/s41388-019-1023-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, J. S., Tachibana, I., Passe, S. M., Huntley, B. K., Borell, T. J., Iturria, N., et al. (2001). PTEN Mutation, EGFR Amplification, and Outcome in Patients with Anaplastic Astrocytoma and Glioblastoma Multiforme. JNCI J. Natl. Cancer Inst. 93 (16), 1246–1256. doi:10.1093/jnci/93.16.1246

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, A., Ra, S., Li, X., Zhou, J., and Binder, S. (2013). Differentiating Cutaneous Squamous Cell Carcinoma and Pseudoepitheliomatous Hyperplasia by Multiplex qRT-PCR. Mod. Pathol. 26 (11), 1433–1437. doi:10.1038/modpathol.2013.82

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, A. A., Fisher, J. L., Rahme, G. J., Hampton, T. H., Baron, U., Olek, S., et al. (2015). Regulatory T Cells Are Not a Strong Predictor of Survival for Patients with Glioblastoma. Neuro-Oncology 17 (6), 801–809. doi:10.1093/neuonc/nou363

PubMed Abstract | CrossRef Full Text | Google Scholar

Töpfer, N., Fuchs, L. M., and Aharoni, A. (2017). The PhytoClust Tool for Metabolic Gene Clusters Discovery in Plant Genomes. Nucleic Acids Res. 45 (12), 7049–7063. doi:10.1093/nar/gkx404

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinuesa, C. G., Linterman, M. A., Yu, D., and MacLennan, I. C. M. (2016). Follicular Helper T Cells. Annu. Rev. Immunol. 34, 335–368. doi:10.1146/annurev-immunol-041015-055605

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Green, M., Choi, J. E., Gijón, M., Kennedy, P. D., Johnson, J. K., et al. (2019). CD8+ T Cells Regulate Tumour Ferroptosis during Cancer Immunotherapy. Nature 569 (7755), 270–274. doi:10.1038/s41586-019-1170-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Wang, Y., Yang, T., Xing, H., Wang, Y., Gao, L., et al. (2021). 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. Bioinform 22 (5), bbab032. doi:10.1093/bib/bbab032

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, P., Pattarini, R., Rong, Y., Guo, H., Bansal, P. K., Kusnoor, S. V., et al. (2012). The Cbln Family of Proteins Interact with Multiple Signaling Pathways. J. Neurochem. 121 (5), 717–729. doi:10.1111/j.1471-4159.2012.07648.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Minikes, A. M., Gao, M., Bian, H., Li, Y., Stockwell, B. R., et al. (2019). Intercellular Interaction Dictates Cancer Cell Ferroptosis via NF2-YAP Signalling. Nature 572 (7769), 402–406. doi:10.1038/s41586-019-1426-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, D., Zhou, Y., Wang, X., Zhao, H., Nie, C., and Jiang, X. (2021). A Ferroptosis-Related Prognostic Risk Score Model to Predict Clinical Significance and Immunogenic Characteristics in Glioblastoma Multiforme. Oxidative Med. Cell. Longev. 2021, 9107857–9107930. doi:10.1155/2021/9107857

CrossRef Full Text | Google Scholar

Xiong, Y., Xiong, Z., Cao, H., Li, C., Wanggou, S., and Li, X. (2020). Multi-dimensional Omics Characterization in Glioblastoma Identifies the Purity-Associated Pattern and Prognostic Gene Signatures. Cancer Cell Int. 20, 37. doi:10.1186/s12935-020-1116-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, H., Parsons, D. W., Jin, G., McLendon, R., Rasheed, B. A., Yuan, W., et al. (2009). IDH1andIDH2Mutations in Gliomas. N. Engl. J. Med. 360 (8), 765–773. doi:10.1056/NEJMoa0808710

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, F., He, Z., Duan, H., Zhang, D., Li, J., Yang, H., et al. (2021). Synergistic Immunotherapy of Glioblastoma by Dual Targeting of IL-6 and CD40. Nat. Commun. 12 (1), 3424. doi:10.1038/s41467-021-23832-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W.-H., Ding, C.-K. C., Sun, T., Rupprecht, G., Lin, C.-C., Hsu, D., et al. (2019). The Hippo Pathway Effector TAZ Regulates Ferroptosis in Renal Cell Carcinoma. Cell Rep. 28 (10), 2501–2508. doi:10.1016/j.celrep.2019.07.107

PubMed Abstract | CrossRef Full Text | Google Scholar

Yee, P. P., Wei, Y., Kim, S.-Y., Lu, T., Chih, S. Y., Lawson, C., et al. (2020). Neutrophil-induced Ferroptosis Promotes Tumor Necrosis in Glioblastoma Progression. Nat. Commun. 11 (1), 5424. doi:10.1038/s41467-020-19193-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Du, M., and Lu, L. (2022). A Novel 16-Genes Signature Scoring System as Prognostic Model to Evaluate Survival Risk in Patients with Glioblastoma. Biomedicines 10 (2), 317. doi:10.3390/biomedicines10020317

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, L., Zhou, X., Xie, F., Zhang, L., Yan, H., Huang, J., et al. (2022). Ferroptosis in Cancer and Cancer Immunotherapy. Cancer Commun. 42 (2), 88–116. doi:10.1002/cac2.12250

CrossRef Full Text | Google Scholar

Zheng, H., Song, K., Fu, Y., You, T., Yang, J., Guo, W., et al. (2021). An Absolute Human Stemness Index Associated with Oncogenic Dedifferentiation. Brief. Bioinform 22 (2), 2151–2160. doi:10.1093/bib/bbz174

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, B., Lawrence, T., and Liang, Y. (2021). The Role of Plasmacytoid Dendritic Cells in Cancers. Front. Immunol. 12, 749190. doi:10.3389/fimmu.2021.749190

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuo, S., Chen, Z., Yang, Y., Zhang, J., Tang, J., and Yang, K. (2020). Clinical and Biological Significances of a Ferroptosis-Related Gene Signature in Glioma. Front. Oncol. 10, 590861. doi:10.3389/fonc.2020.590861

PubMed Abstract | CrossRef Full Text | Google Scholar