Validated ferroptosis driver and suppressor genes were collected from the FerrDb (Table S1, http://www.zhounan.org/ferrdb/legacy/operations/download.html) (8). The TCGA RNA-Seq data of primary tumor samples from 21 cancer types were downloaded from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). Overall survival (OS), age, gender, and race information of each TCGA cancer was obtained from the integrated TCGA Pan-cancer clinical data resource (9). Abbreviations and patient numbers of each TCGA cancer type are included in Table S2. Independent datasets GSE104462, GSE121689, GSE32369, GSE132257, GSE103322, GSE123904, GSE87211, GSE83129, GSE34228, GSE16648, GSE81005, and GSE156632 were available on GEO (http://www.ncbi.nlm.nih.gov/geo , Table S3). The Clinical Proteomic Tumor Analysis Consortium (CPTAC) CRC and LUAD protein expression data were obtained from https://proteomics.cancer.gov/programs/cptac.

According to instruction from FerrDb (8), drivers, which are genes that promote ferroptosis, and suppressors, which are genes that prevent ferroptosis. Then, ferroptosis driver genes and suppressor genes were used to calculate drive and suppress the activity by ssGSEA, respectively. Thus, drive activity represents the level of promoting the ferroptosis biological process and suppress activity represents the level of preventing the ferroptosis biological process. Finally, we defined   R F L = Drive activity Suppress activity We estimated RFL-correlated genes in the tumors of nine cancers from the TCGA. A total of 19,211 protein-coding genes were selected in the analysis. Those with a correlation between gene expression and RFL of R >|0.3| and P <0.05 were considered markedly correlated. The correlation levels were based on the Pearson test. According to a previous study (10), the Jaccard index was used to evaluate the ratio of the RFL-correlated genes common to the two cancer types (10).

CRC (RKO, HCT116, and HCT8) and LUAD (A549 and H1299) cell lines were cultured as described in previous studies (11, 12). Erastin (MedChemExpress, China) was used to induce ferroptosis in these cancer cell lines. Total RNA was extracted using TRIzol reagent (Takara, Japan). Complementary cDNA was generated using the Prime-Script RT reagent kit (CWBIO, China). The UltraSYBR Mixture (CWBIO, China) was used to detect the relative mRNA expression. The levels of β-actin were used as the reference and normalized control. The primers are listed in Table S4.

Four different methods were used to evaluate immune infiltrates in TCGA tissue-level expression data, namely, ESTIMATE, XCELL, 22 immune cell types (LM22) of CIBERSORTx, and Pan-cancer immunogenomic based on single sample gene set enrichment analysis (ssGSEA) (13–18). The Seurat package in the R software was used to perform single cell cluster and annotation (19). Cells from each patient were annotated as related studies of GSE132257, GSE103322, GSE123904, and GSE156632 (20–22). The Scissor method was used to identify ferroptotic status at single-cell level based on RFL in each TCGA cancer type (23). A TIDE score was calculated online (http://tide.dfci.harvard.edu/) to compute for tumor samples to predict patient response.

OS, progression-free survival (PFS), and relapse-free survival (RFS) were evaluated by Kaplan–Meier survival analysis and log-rank test as previously described (10, 18). For Figures 2A–C, the clinical samples of each cancer type we used were divided into two groups with high and low RFL group based on the level of RFL. All RFL values from the 10th to 90th percentiles were used as cutoffs to group the samples for survival analyses, respectively; the value yielding the lowest log-rank P -value was selected as the cutoff value to group samples. For Figure 3B, the clinical samples of each cancer type we used were divided into two groups with high and low expression based on the expression level of each select gene, respectively. All expression values from the 10th to 90th percentiles were used as cutoffs to group the samples for survival analyses, respectively; the value yielding the lowest log-rank P -value was selected as the cutoff value to group samples. Hazard ratio (HR) >1 and P < 0.05 were considered significant and associated with poor survival. Hazard ratio (HR) <1 and P <0.05 were considered significant and associated with favorable survival. The survival of patients was analyzed by R software 4.1.0 and Graphpad prism 9. Cox regression model analyses were conducted using Graphpad prism 9. T test or Mann–Whitney test was used in the two-group comparison, and the Kruskal–Wallis test was used in the three-group comparison using Graphpad Prism 9. All reported P-values were two-sided. Group comparison, receiver operating characteristics (ROC) curve testing, and Chi-squared test analyses were performed using GraphPad Prism 9. Analyses were performed using R software 4.1.0 with the GSVA (1.40.1) package (ssGSEA) for Pan-cancer immunogenomic and enrichment score calculation; survival (3.2–11) package for two-group OS analysis; ESTIMATE (1.0.13) package for ESTIMATE ImmuneScore and TumorPurity estimation; Seurat (4.0.3) package for cluster and annotation in single-cell GC samples; and Scissor (2.0.0) package for identifying of ferroptotic status at single-cell level by integrating bulk sequencing data. A P-value of less than 0.05 was regarded as statistically significant. Figures were designed, analyzed, and visualized by GraphPad Prism 9 and R software 4.1.0.

To classify the ferroptosis level of tumor samples, we focused on the FerrDb, the first database of ferroptotic regulator collection in the world (Figure 1A) (8). First, validated ferroptosis driver genes (n = 86) and suppressor genes (n = 66) were collected according to the FerrDb. Four genes were excluded because they belong to both the ferroptosis driver and suppressor gene signatures. Next, protein-coding genes were selected for further mutation analysis for reliability. The mutation rates of these genes were lower than 5%, except for PIK3CA (13%). Then, 79 ferroptosis driver genes and 55 suppressor genes were used to calculate the drive and suppress the activity by ssGSEA, respectively. Finally, RFL was established based on drive active and suppressed active.

Multiple analyses were performed to validate the performance and robustness of this RFL. We used the TCGA pan-cancer data to establish RFL and found positive correlations between RFL and drive active and negative correlations between RFL and suppress active in 21 cancer types (Figure 1B). Moreover, our data in the GEO datasets showed that cancer cells under ferroptosis inducer Erastin conditions showed significantly higher RFL, whereas those under ferroptosis inhibitor DFO conditions showed significantly lower RFL (Figure 1C). These results demonstrate the robustness of the RFL to define and evaluate ferroptosis levels in human cancers.

To assess the relevance of ferroptosis level in cancer clinical trials, we examined the correlations of our RFL classification with the survival time of patients in the TCGA human cancers according to our previous pan-cancer studies (10). LUSC, PAAD, and THCA were excluded with correlations between RFL and drive or suppress active <|0.2|. Interestingly, high FRI predicts favorable OS and PFS in most cancer types (Figure 2A). We further concentrated on OS and found that high- and low-RFL groups greatly varied across different cancer types based on OS (Figure 2B). Moreover, we observed that high RFL cases were consistently associated with a significantly favorable prognosis across CESC, CHOL, CRC, GBM, HNSC, KIRC, KIRP, LGG, and LUAD (Figure 2C). Importantly, potential confounding factors, such as age, gender, race, and TNM stage, were analyzed with RFL to evaluate the independent effect of COX models. High RFL remains a significant independent factor that predicts favorable survival in CHOL, CRC, GBM, HNSC, KIRC, KIRP, LGG, and LUAD, though P = 0.0557 in CESC (Supplementary Figure 1). Additionally, the high RFL group demonstrated higher ferroptosis drive active in these cancer types, whereas the low RFL group demonstrated lower ferroptosis suppress active (Figure 2D). Hence, these results suggest the potential prognostic power of RFL classification, which considers ferroptosis drive and suppress regulators.

To identify new RFL-correlated genes in human cancer, we analyzed correlations between RFL and protein code genes across CESC, CHOL, CRC, GBM, HNSC, KIRC, KIRP, LGG, and LUAD. A series of RFL-positive and -negative genes have been identified in these cancer types, and we used Jaccard indexes to evaluate the overlap of these RFL-correlated genes in each cancer type. RFL-positive genes showed different patterns than RFL-negative genes (Figure 3A). Next, we found 21 RFL-correlated genes were positively (n = 16) or negatively (n = 5) correlated with RFL in at least seven cancer types (Figure 3B). Interestingly, most of these RFL-positive genes exhibit favorable OS in human cancers, while RFL-negative genes show poor survival. These genes were further selected to develop Cox regression models for each cancer type (Figure 3C). The risk scores decreased as the RFL and RFL-positive genes were expressied, while they increased as RFL-negative genes were expressed in most cancer types, except LGG. The high-risk group (risk score of >0) had more death status compared with the low-risk group (risk score of <0). Moreover, our analysis showed that there was a significant difference in prognosis and poor survival between the high and low-risk groups in CESC, CHOL, CRC, GBM, HNSC, KIRC, KIRP, LGG, and LUAD (Figures 3D, E). These data hint that the models according to the RFL-correlated genes identified by RFL are effective in the prognosis of human cancers.

To determine the functional significance of the RFL-correlated genes identified by RFL above, seven candidate genes that have not been reported to be ferroptosis-related were selected and experimentally validated in five cancer cell lines. Most of these genes were upregulated after Erastin-induced ferroptosis in RKO, HCT116, HCT8, A549, and H1299 cell lines (Figure 4A ). Intriguingly, most of these RFL-positive genes were downregulated in tumor tissues at mRNA and protein levels of CRC and LUAD (Figure 4B). Thus, RFL-correlated genes identified by RFL are associated with Erastin-induced ferroptosis and tumorigenesis. Our findings suggest the potential of the RFL, which can be an effective resource to functionally identify RFL-correlated genes.

The above mentioned data demonstrate robustness, potential prognostic, and identify the power of the RFL. We therefore explore the molecular characteristics of different RFL groups based on previous subtype studies performed by the TCGA Research Network, etc. The RFL high group comprised more MSI/CpG island methylator phenotype (CIMP) tumors in CRC (Supplementary Figure 2A). The MSI is mostly located in the right colon and is frequently associated with the CIMP (24). Interestingly, current clinical dogma considers MSI tumors to carry a good prognosis and a low risk of relapse (25). Previous research has exhibited a sustained clinical response to immune checkpoints (ICs) with dramatic clinical improvement in CRC patients with MSI-H (26). In HNSC, the RFL high group showed fewer classical cases compared with the RFL low group (Supplementary Figure 2B). TP53 mutation, alteration of oxidative stress genes, and heavy smoking history occurred in most classical tumors (27). RFL high group showed more proximal-inflammatory subtype and lower proximal-proliferative subtype in LUAD (Supplementary Figure 2C) (28). In KIRC, high RFL cases comprised more m1 and m3 subtypes, while RFL low cases had more m2 and m4 subtypes (Supplementary Figure 2D) (29). High RFL KIRP patients showed a few C2c subtype cases compared with the low RFL group (Supplementary Figure 2E). C2c tumors consisting solely of CIMP-associated KIRP patients, who were younger at the time of presentation and had a lower probability of OS, were also characterized by a Warburg-like metabolic shift to glycolysis-dependent metabolism and upregulation of hypoxia-related genes (30). In LGG, the low RFL group showed a higher proportion of IDH wild-type patients (Supplementary Figure 2F). The large majority of LGG without an IDH mutation had genomic aberrations and clinical behavior strikingly similar to those found in primary GBM (31). In GBM, high RFL cases largely overlap with the mesenchymal subtype, which showed high levels of genes in the tumor necrosis factor superfamily pathway and NF-κB pathway, suggesting that higher overall necrosis and associated inflammatory infiltrate (Supplementary Figure 2G) (32, 33). These data provide an overview of the molecular differences associated with RFL across human cancers, particularly parts of characterization of these molecular subtypes that are associated with RFL level and matched clinical outcomes.

Analyses of the molecular characterization of ferroptosis showed that high RFL was associated with inflammatory infiltrates. Recent studies have suggested correlations between ferroptosis-correlated factors and proinflammatory TME, which activates anti-tumor immune responses, suggesting a unique relationship between ferroptosis and TME (34, 35). To further investigate this, four independent methods were employed to evaluate immune infiltrates in CESC, CHOL, CRC, GBM, HNSC, KIRC, KIRP, LGG, and LUAD (13–17). Notably, RFL was positively correlated with total immune infiltrates and negatively correlated with tumor purity in these cancer types from TCGA, except LGG (Figure 5A).

The availability of single-cell transcriptomic profiles across a broad range of human cancers provides an unprecedented opportunity to explore TME in greater depth. Thus, we used single-cell datasets from GEO to gain a better understanding of ferroptosis-correlated TME. First, fresh cell preparation single-cell data of CRC were used to analyze and we identified cancer cells and a series of nonmalignant cells (Figure 5B). New biotechnology that identifies cell subpopulations from single-cell data according to a given phenotype-associated bulk expression data,was used to identify high- and low-RFL cell subpopulations in this CRC single-cell data (Figure 5C) (23). Interestingly, the high-RFL cell subpopulation contained more nonmalignant cell types (i.e., T cells and macrophages) and fewer cancer cells compared with the low-RFL cell subpopulation (Figure 5D). We also performed similar analyses on HNSC single-cell data and observed more significant results, suggesting that high-RFL has a comparatively fertile ground of immune cell infiltration in TME (Figures 5E, G). Next, we used another method to analyze LUAD single-cell data, and patients with primary tumors were selected. The mean values of single-cell expression in each patient sample for each gene were calculated as pseudobulk transcriptome data to obtain RFL at the tissue level of each case (Figure 5H). The LX661 and LX680 patients showed the highest and lowest RFL, respectively (Figure 5I). Then, we identified cancer cells and a series of nonmalignant cells in these two samples. Impressively, the high-RFL (LX661) sample exhibited more T-cell infiltration and lower tumor purity compared with the low-RFL (LX680) sample (Figure 5J). Furthermore, RFL was downregulated in cancer cells compared with nonmalignant cell types (Figure 5K). We also analyzed the KIRC single-cell data of primary tumors and observed similar results (Figures 5L–O). These results hint that RFL reflects immune cell infiltration and tumor purity in the TME of human cancers. The level of ferroptosis should be considered and could be a potential strategy to provide immunostimulatory in human cancers. This strategy has been reported in a new study of nanomedicine, which demonstrates simultaneous targeting of both the TME and cancer cells (36).

To characterize the clinically applicable therapeutic implications of RFL, we first examined cancer immunotherapy due to the abovementioned data of TME. Immunomodulators play a critical role in cancer immunotherapy and we found that most of the immunomodulator genes were overexpressed in high RFL tumors, except LGG (Supplementary Figure 3) (37). Next, nine ICs were selected according to the previous TME study, and RFL was positively correlated with these IC expressions in CESC, CHOL, CRC, GBM, HNSC, KIRC, KIRP, and LUAD from the TCGA, except LGG (Figure 6A) (38). These results suggest that high RFL is correlated with an immunosuppressive phenotype. Besides, RFL was negatively correlated with the TIDE score, which is computed for tumor samples and can serve as a surrogate biomarker to predict the response to ICs blockade (Figure 6B). The patient who have low TIDE might represent the response to IC inhibitor-based treatment (39). However, CRC showed a weak positive correlation between RFL and TIDE. These data on TME, ICs and response prediction motivated our next analysis in immune overdrive subtype, a novel TME-based immune molecular subtyping system have reported in previous research (18). The FOXP3^HighCTLA4^High* subpopulation (Treg subpopulation with overexpressed CTLA4) was characterized by an immune overdrive phenotype with high immune cell infiltration, low tumor purity, high IC levels, and TGF-β activation in CRC, and observed in LGG. The FOXP3^HighCTLA4^High* cases showed higher RFL compared with other groups in the CRC (Figure 6C). Nevertheless, we noted that the FOXP3^HighCTLA4^High* cases have the lowest RFL compared with other groups in LGG, consistent with the abovementioned data. These results suggest that patients with tumors have a relatively higher ferroptosis level, which may be related to combination treatment with immunotherapy, particularly IC inhibitor combination treatment.

Traditional chemotherapy or radiotherapy remains the selection of first-line therapy in multiple cancer types. Thus, we first evaluated the effect of ferroptosis status in preoperative radiochemotherapy of CRC. High RFL was associated with a lower risk of recurrence in CRC patients with no matter 5-FU + radiotherapy or 5-FU + oxaliplatin + radiotherapy (Figure 6D). Responder cases treated first-line with 5-FU and oxaliplatin demonstrated higher RFL than non-responder cases (Figure 6E). Chemotherapeutic drug-sensitive cancer cell lines, such as gefitinib (PC9), methotrexate (HT29) or 5-FU (HCT8), show a higher RFL compared with matched drug-resistant cell lines (Figure 6F). Taken together, our analyses show that the ferroptosis status meaningfully correlated with tumor response to immunotherapy, chemotherapy, or radiotherapy. Thus, the tumor ferroptosis status should be considered to improve the efficacy of these cancer therapies.