Development and Validation of a Combined Ferroptosis and Immune Prognostic Classifier for Hepatocellular Carcinoma

Background Immunotherapy and sorafenib exert anti-tumor effects via ferroptosis, but reliable biomarkers for the individual treatment and prognosis prediction of hepatocellular carcinoma (HCC) based on the ferroptosis and immune status remain lacking. Methods Ferroptosis-related genes (FRGs) were identified by downloading data from FerrDb and by searching and reading original articles from PubMed. Immune-related genes (IRGs) were downloaded from ImmPort. Prognostic FRGs and IRGs in the GSE14520 (n = 220) and The Cancer Genome Atlas (TCGA, n = 365) datasets were identified. Least absolute shrinkage and selection operator (LASSO) Cox regression and multivariate Cox regression were used for model construction. Ferroptosis expression profiles, the infiltration of immune cells, and the somatic mutation status were analyzed and compared. Results Twenty-seven prognostic ferroptosis- and immune-related signatures were included to construct a comprehensive index of ferroptosis and immune status (CIFI). A subgroup of patients was identified as having a high CIFI value, which was associated with a worse prognosis. This subgroup of patients had significantly up-regulated expressions of many suppressors of ferroptosis and higher fractions of immunosuppressive cells, such as cancer-associated fibroblasts (CAFs) and myeloid-derived suppressor cells (MDSCs). Notably, somatic mutation analysis indicated that high-CIFI patients had higher levels of tumor heterogeneity and higher mutation frequencies of genes like TP53. Conclusion In this work, a novel prognostic classifier was developed based on ferroptosis- and IRGs in HCC, and this classifier could be used for prognostic prediction and the selection of patients for immunotherapies and targeted therapies.


INTRODUCTION
Hepatocellular carcinoma (HCC) is a lethal disease with high incidence and dismal prognosis (Hartke et al., 2017;Kulik and El-Serag, 2019). Even among patients who receive resection at an early stage, HCC reoccurs in more than half within 2 years (Tabrizian et al., 2015). Sorafenib and lenvatinib are the two first-line therapies for advanced HCC patients, who show a medium overall survival (OS) of around 1 year with the treatment (Llovet et al., 2008;Kudo et al., 2018). Many phase I and II clinical trials suggest that immunebased therapies, including anti-CTLA-4, anti-PD-1, and anti-PD-L1 strategies, benefit HCC patients (Hilmi et al., 2019). A recent phase III study also indicated that unresectable HCC patients who receive atezolizumab plus bevacizumab as a firstline treatment have a 1-year survival rate of 67.2% and a median progression-free survival of 6.8 months, which are better than those of patients who receive only sorafenib (Finn et al., 2020a). However, another immune checkpoint inhibitor, pembrolizumab, failed to reach its primary end-points as a second-line therapy against HCC in a randomized, double-blind, phase III trial (Finn et al., 2020b). Nivolumab also failed to prove its superiority over sorafenib in a phase III trial . These results suggest that HCC might have a complex immune status, and more studies are required to understand its underlying mechanisms.
Ferroptosis might be one such essential mechanism in HCC that deserves attention. As an iron-dependent, lipid peroxidation-mediated form of cell death, ferroptosis is distinct from apoptosis, necrosis, and autophagy (Dixon et al., 2012). Since its first definition in 2012, numerous genes have been identified to regulate this new form of cell death (Dixon et al., 2012;Hassannia et al., 2019). Some are well-defined suppressors of ferroptosis (SOFs), some are drivers of ferroptosis (DOFs), and some genes, like TP53, could act as a SOF or DOF in a context-dependent manner (Jiang et al., 2015;Tarangelo et al., 2018). A remarkable study was conducted by Wang W. et al. (2019), who reported that CD8+ T cells enhance ferroptosis by down-regulating SLC3A2 and SLC7A11, and the induction of ferroptosis contributes to the anti-tumor efficacy of immunotherapy, suggesting that the immune system might, at least in part, function through ferroptosis (Stockwell and Jiang, 2019). Moreover, ferroptotic cancer cells might release signals like oxidized lipid mediators to affect anti-tumor immunity, or a small proportion of cells undergoing ferroptosis might suppress the immune system and allow tumor growth (Friedmann Angeli et al., 2019). However, a comprehensive regulatory network between ferroptosis and immune response is currently lacking, as few studies have explored their relationship.
Consequently, ferroptosis, a biological process in which immune therapy and sorafenib converge to exert an anti-tumor effect, might have a fundamental impact on the treatment response of HCC. The understanding of the occurrence and regulation of ferroptosis has increased considerably in recent decades; however, only a few studies have explored the ferroptosis-related genes (FRGs) and pathways in HCC (Louandre et al., 2013(Louandre et al., , 2015Carlson et al., 2016;Houessinon et al., 2016;Sun et al., 2016b;Bai et al., 2019;Qi et al., 2019;Feng et al., 2020;Kim et al., 2020). With the currently available FRGs and immune-related genes (IRGs), and the accumulative data deposited in public databases like The Cancer Genome Atlas (TCGA), it is hypothesized that a prognostic molecular classifier based on the immune response and ferroptosis status might help to identify subgroups of HCC patients with distinct ferroptosisimmune phenotypes and survival profiles. In the present work, it is demonstrated that a comprehensive index of ferroptosis and immune status (CIFI) developed from FRGs and IRGs is tightly correlated with the actual ferroptosis and immune status, as well as the prognosis of HCC patients.

Study Population and Data Acquisition
Publicly available data on HCC were downloaded from the Gene Expression Omnibus (GEO) database (GSE14520/GPL3921) and TCGA, and were processed as reported in the authors' previous study . Patients who met the following selection criteria were included: (a) histologically diagnosed with HCC; (b) available gene expression data; (c) available survival information. The first sample according to the label was selected if the same patient had two or more samples in the TCGA dataset. The baseline characteristics of patients in these two cohorts are summarized in Table 1.

Generation of FRGs and IRGs
The FRGs were first downloaded from the FerrDb website. 1 Because the database was last updated on February 23, 2020, PubMed 2 was searched with the keyword "ferroptosis" on July 13, 2020, and 1,113 records were obtained. The authors read all the original articles and supplemented the list of FRGs (Supplementary Table 1) with the following criteria: (a) the inhibitor(s) or activator(s) of the gene could regulate ferroptosis; (b) the over-expression or down-regulation of the gene could regulate ferroptosis. Non-coding RNAs (ncRNAs) and markers of ferroptosis were not included. A total of 283 FRGs were identified, and are summarized in Supplementary Table 1. These FRGs were further divided into two groups based on their regulatory effects on ferroptosis, namely DOFs and SOFs.  The comprehensive list of IRGs was downloaded from ImmPort 3 on August 2, 2020, as reported in other studies Zheng et al., 2020). The FRGs and IRGs available in the GSE14520 and TCGA datasets were included in this study.

Development of the CIFI
The FRGs and IRGs in the GSE14520 and TCGA were first respectively subjected to univariate Cox regression. The genes with prognostic significance (p < 0.05) in both datasets were input into a least absolute shrinkage and selection operator (LASSO) Cox regression model, which is a widely used machine learning algorithm that deals with multicollinearity (Gui and Li, 2005;Wu et al., 2019). The analysis generated crucial genes for model construction, and was achieved by using the "glmnet" package in R software (Friedman et al., 2010). A multivariate Cox regression model was applied to obtain the regression coefficients for these crucial genes. To reflect the comprehensive effects of the ferroptosis and immune status, a risk score of each patient was calculated by multiplying the normalized gene expression of each crucial gene with its corresponding regression coefficient. To facilitate the interpretation of results from different datasets, the CIFI value of each patient was calculated by using the risk score of the patient subtracted by the minimum risk score of the cohort, which was then divided by the maximum risk score of the cohort, namely CIFI = (Risk score−Min)/Max.

Functional Analysis and Heterogeneity
Gene Set Enrichment Analysis (GSEA software, version 4.0.3) was used to investigate the pathways enriched in the high-and low-CIFI subgroups. Gene expression data were loaded into GSEA, and c2.cp.v7.1.symbols.gmt was selected as the gene set database. The pathways with the following criteria were regarded to be significantly enriched: nominal p-value < 0.05, false discovery rate (FDR) q-value < 0.25, and normalized enrichment score (NES) > 1.
The crucial genes were constructed into a protein-protein interaction (PPI) network by uploading them into the STRING database. 4 The network was visualized by Cytoscape (version 3.8.0). The expressions of these genes in tumor and normal samples were visualized by the "pheatmap" package in R. Principal component analysis (PCA) was carried out to examine the clustering efficacy of the prognostic signature with the "gmodels" package in R.
Somatic mutation data of the TCGA cohort were downloaded from the GDC database on August 7, 2020. The downloaded MAF files of simple nucleotide variation (workflow type: varScan2 variant aggregation and masking) were processed and visualized by the "maftools" package in R. The tumor mutation burden (TMB) and the mutant-allele tumor heterogeneity (MATH) score of tumor samples in the TCGA cohort were also calculated by the "maftools" package.

Immune Profile Analysis
To analyze the immune status of each sample in the GSE14520 and TCGA cohorts, the relative infiltrations of 28 immune cell types in the tumor microenvironment (TME) were calculated via single-sample GSEA (ssGSEA) with the application of the "GSVA" package in R. The feature gene panels for each type of immune cell were downloaded from the publication by Charoentong et al. (2017). The normalized gene expression data of the GSE14520 and TCGA cohorts were further uploaded into the EPIC website, 5 and the proportions of eight categories of immune cells were estimated according to the instructions on the website (Racle et al., 2017).

Statistical Analysis
Univariate and multivariate Cox regressions were conducted by using the "survminer" package in R. The OS, relapse-free survival (RFS), disease-specific survival (DSS), and progressionfree interval (PFI) of the high-and low-CIFI subgroups were compared using the Kaplan-Meier method with a log-rank test. Time-dependent receiver operating characteristic (ROC) analyses and the subsequent calculation of the area under the curve (AUC) were performed using the "timeROC" package in R. The "corrr" and "corrplot" packages were used to conduct the correlation analysis with the Pearson method. The significance of the differences in the expressions of specific genes or in the fractions of immune cells was assessed by the Wilcoxon test. Student's t-test was used to compare the differences in the TMB and MATH scores between the highand low-CIFI subgroups. The "ggplot2, " "ggforest, " "cowplot, " "plot3D, " "maftools, " "VennDiagram, " and "ggplotify" packages in R (version 4.0.2) were used for visualization. A p-value of less than 0.05 was considered to be statistically significant ( * p < 0.05; * * p < 0.01; * * * p < 0.001; * * * * p < 0.0001).

Construction of the CIFI Classifier in HCC
Populations of 220 patients from the GSE14520 dataset and 365 patients from the TCGA dataset were identified and included in this study ( Table 1). The gene expression data were used to construct the CIFI classifier. In total, 279 FRGs (four genes were removed because they were identified as both DOFs and SOFs) and 952 IRGs in the two datasets were mapped and used for model construction ( Figure 1A). Univariate Cox regression analysis was conducted to estimate the prognostic significance of these genes. The results showed that 85 genes in both the GSE14520 and TCGA cohorts had significant prognostic relevance (p < 0.05) ( Figure 1B and Supplementary Table 2). The expressions of these 85 genes in the GSE14520 cohort are shown in Figure 1C. The PPI network of these genes demonstrates that the FRGs had strong correlations with the IRGs (Figure 1D).
The LASSO Cox regression model was then applied to construct a prognostic model for the OS of HCC patients by using the gene expression data of the 85 genes in the GSE14520 dataset. The model identified 27 genes based on the optimal value of λ (Figures 2A,B). Among them, G6PD, WDR76, CA9, and AHCYL1 were FRGs, HMOX1, and FLT3 participated in both ferroptosis and the immune process, and the remaining 21 genes (SPP1, EPO, CKLF, GLP1R, LHB, NR1H3, ADM, GAL, IRF5, IL18RAP, SEMA3F, PLXNA2, MMP12, ANGPT1, ECD, FABP3, LANCL1, OGFR, GH2, STC1, and OSMR) were IRGs. A novel risk score was calculated by multiplying the gene expression of each gene and its corresponding coefficient, which was obtained by multivariate Cox regression analysis. The CIFI values were generated by the formula mentioned in section "Materials and Methods." The correlation analysis indicated that the CIFI value was significantly correlated with the selected genes ( Figure 2C), and the correlations between the CIFI values and FRGs are shown in Figure 2D. The HCC patients were stratified into high-risk (high-CIFI, n = 80) and low-risk (low-CIFI, n = 140) subgroups based on the optimal CIFI cut-off value (0.52), which was calculated by the surv_cutpoint function in the "survminer" package ( Figure 2E and Supplementary Figure 1A). As illustrated in Figure 2F, the patients in the high-risk subgroup  The LASSO Cox regression model was constructed from the 85 prognostic genes, and the tuning parameter (λ) was calculated based on the partial likelihood deviance with 10-fold cross-validation. An optimal log λ value is indicated by the vertical black line in the plot. (C,D) Correlation networks (C) between the CIFI value and the 27 signature genes or (D) between the CIFI value and the ferroptosis-related genes in the GSE14520 dataset. (E,F) The distribution and optimal cutoff value of (E) the risk scores and (F) the OS status and OS in the GSE14520 dataset. (G,H) The (G) 2D and (H) 3D plots of the PCA of the GSE14520 dataset based on the expression profiles of the 27 signature genes in different risk groups.
Frontiers in Cell and Developmental Biology | www.frontiersin.org had more occurrences of death and shorter survival times. PCA revealed that the patients in these two subgroups were distinctively clustered (Figures 2G,H).
Time-dependent ROC curves were plotted by R software, and the AUC was calculated at different time points to estimate the predictive performance of the CIFI. As shown in Figure 3A, the AUC reached 0.76 at 1 year, 0.72 at 3 years, and 0.77 at 5 years, suggesting a favorable predictive value of the CIFI in shortand long-term follow-up. Kaplan-Meier curves indicated that patients with high CIFI values had significantly shorter OSs than their counterparts with low CIFI values ( Figure 3B, p < 0.001). In addition, patients with high CIFI values relapsed earlier than those with low CIFI values ( Figure 3C, p < 0.001).

Validation of the CIFI Classifier in HCC
To validate the indicative CIFI value in a larger cohort of HCC patients, the CIFI values in the TCGA dataset (n = 365) were calculated using the same risk formula and cutoff point obtained from the GSE14520 dataset. Ultimately, 63.64% of the HCC patients (n = 232) in the TCGA cohort were categorized into the low-risk (low-CIFI) subgroup, while the remaining patients (n = 133) were categorized into the high-risk (high-CIFI) subgroup. Consistent with the results of the GSE14520 dataset, PCA indicated that the two subgroups in the TCGA cohort were distributed in discrete directions (Figures 4A,B). The AUCs for OS were 0.69 at 1 year, 0.7 at 3 years, and 0.73 at 5 years, indicating an increasing predictive value in long-term follow-up ( Figure 4C). The medium OS time of the patients in the low-risk subgroup was 2,532 days, which was dramatically longer than that of patients in the high-risk subgroup (medium OS time: 899 days, p < 0.0001, Figure 4D). In addition, patients with low CIFIs had significantly longer DSSs (p = 0.00024, Figure 4E) and PFIs (p = 0.0038, Figure 4F).

Ferroptosis Profile in the CIFI
A correlation analysis of the GSE14520 dataset revealed a correlation between the CIFI values and FRGs ( Figure 2D).
Because genes might either facilitate or suppress ferroptosis, the transcriptional changes of SOFs and DOFs in the CIFI-stratified subgroups were investigated. ATF4, CA9, EGLN1, ELAVL1, FTH1, GPX4, HELLS, ITGB8, NFE2L2, SLC7A11, SQSTM1, and VDAC2 are well-investigated SOFs (Dixon et al., 2012;Sun et al., 2016b;Chen et al., 2017;Jiang et al., 2017;Zhang et al., 2018;Bai et al., 2019;. It was found that, excluding GPX4, NFE2L2, and SQSTM1, the remaining SOFs were significantly up-regulated in the high-CIFI subgroup of the GSE14520 cohort ( Figure 5A). The change of these genes was further validated in the TCGA cohort ( Figure 5B). Although GPX4 was found to be significantly down-regulated in the high-CIFI subgroup of the GSE14520 dataset, such a change was not observed in the TCGA cohort (Figures 5A,B). Moreover, it was found that there were no significant differences in most of the DOFs (ALOX12, ALOX12B, ALOX15, ALOXE3, BECN1, BID, GOT1, PRKAA1, and PTGS2) at the transcriptional level between the high-and low-CIFI subgroups in the GSE14520 cohort (Supplementary Figure 1B). Moreover, ALOX15B and NOX1 were significantly downregulated in the high-CIFI subgroup, while ACSL4 and ALOX5 were significantly up-regulated (Supplementary Figure 1B).
Gene Set Enrichment Analysis was further conducted between the two subgroups of both the GSE14520 and TCGA cohorts. The results showed that the genes in the low-CIFI subgroups of the two cohorts were significantly enriched in ferroptosisrelated biological processes like biological oxidation, fatty acid metabolism, peroxisome, mitochondrial fatty acid betaoxidation, peroxisomal lipid metabolism, glyoxylate metabolism, and glycine degradation ( Figure 5C). Considered together, these results suggest a ferroptosis-suppressed status of the high-CIFI subgroup.

Immune Profile in CIFI
To understand how the CIFI reflects the immune status of HCC, ssGSEA was first employed to calculate the immune enrichment scores of various immune categories in each patient, and their relationships with the CIFI value were  investigated. Myeloid-derived suppressor cells (MDSCs) are immune-suppressing cells that contribute to the growth and invasion of HCC, while eosinophil inhibits the growth of the disease (Lin et al., 2017;Hollande et al., 2019;Lu et al., 2019). As shown in Figure 6A, the CIFI value was found to have a significant positive correlation with MDSCs in both the GSE14520 and TCGA datasets, while a significant negative correlation was observed between the CIFI value and eosinophil in both datasets. However, the CIFI value also exhibited positive associations with activated CD4 T cells and central memory CD4 T cells (Figure 6A), which are regarded as cytotoxic cells against HCC (Jin et al., 2018). The proportion of immune cells was also estimated by the EPIC application, and it was found that the CIFI value was positively correlated with the fraction of cancer-associated fibroblasts (CAFs) and negatively correlated with the fraction of macrophages (Figures 6B,C). The GDC website 6 contains data on the fractions of immune cells in tumors of TCGA via application of the CIBERSORT method. As shown in 6 https://gdc.cancer.gov/about-data/publications/panimmune Figure 6D, macrophages M0, plasma cells, neutrophils, and regulatory T cells (Tregs) were up-regulated in the high-CIFI subgroup of the TCGA cohort, while naïve B cells, resting mast cells, monocytes, resting natural killer (NK) cells, naïve CD4 T cells, and CD8 T cells were significantly downregulated (p < 0.05).

Independent Prognostic Value of the CIFI Classifier
The preceding analysis suggests that HCC patients with high CIFI values had a ferroptosis-suppressive and immune-suppressive status; this spurred an interest to analyze the associations between the CIFI value and clinicopathological features of these patients. As shown in Table 2, more patients in the high-CIFI subgroup in both the GSE14520 and TCGA cohorts were at an advanced TNM stage (stage III or IV, p < 0.001), and a considerably higher percentage of patients in the high-CIFI subgroup of the GSE14520 cohort were at an advanced BCLC stage (p = 4.88 × e −5 ). In addition, a significantly higher percentage of HCC patients in the high-CIFI subgroup of the GSE14520 cohort had cirrhosis (p = 0.009881), while no difference was observed in patients with fibrosis between the two subgroups of the TCGA cohort (p = 0.2979). Moreover, more patients in the high-CIFI subgroup of the TCGA cohort had higher histologic grades (G3 or G4, p < 0.001) or vascular invasion (p = 0.02296). High-CIFI HCC patients in the GSE14520 cohort had higher levels of AFP (p = 0.001172), and such a relationship was also observed in the TCGA cohort, although the result did not reach significance. No consistent differences were observed between CIFI value and age or gender in the GSE14520 and TCGA cohorts. Considered together, these results suggest that a high CIFI value is associated with worse clinicopathological features of HCC patients.
To determine whether the CIFI could serve as an independent prognostic predictor of OS, the clinicopathological features and CIFI values were first input into a univariate Cox regression analysis. As illustrated in Figures 8A, 9A, the CIFI value was found to be significantly associated with OS in both the GSE14520 (HR = 7.785, 95% CI = 4.848-12.5, p < 0.001) and TCGA datasets (HR = 2.498, 95% CI = 1.763-3.541, p < 0.001). Then, the CIFI values and the clinicopathological features with prognostic significance (p < 0.05) were subjected to a multivariate Cox regression analysis, which revealed that CIFI remained an independent prognostic predictor after correction for other confounding factors (p < 0.001, Figures 8B, 9B).

DISCUSSION
The advent of selective genome-and immune-targeted therapy has considerably improved the prognosis of cancer patients, and treatment decisions tend to be made on the basis of abnormal molecular profiles or "signatures" rather than the tissue type or anatomical site of tumors (Jackson and Chester, 2015;Abubakar and Gan, 2016;Jia et al., 2016). For instance, patients with lung cancer, colon cancer, or melanoma bearing BRAF V600E mutations all benefit from inhibitors targeting this mutation (Long et al., 2017;Planchard et al., 2017;Kopetz et al., 2019). Indeed, an increasing number of studies have attempted to identify subgroups of tumor patients based on their molecular profiles, which reflect unique phenotypes, distinct treatment responses, and different prognoses. For example, the hypoxia low /immune high subgroup of triple-negative breast cancer (TNBC) patients has a significantly longer OS and might have a better response toward immunotherapies . Colorectal cancer (CRC) patients can be classified into low-and high-risk groups on the basis of their autophagy-related features, and high-risk CRC patients might require more aggressive treatment interventions .
Considering that immunotherapy and sorafenib inhibit the growth of tumors (including HCC) via the induction of ferroptosis (Louandre et al., 2013;Wang W. et al., 2019), HCC patients with distinct ferroptosis and immune phenotypes might have different prognoses. In this work, the CIFI of HCC was developed based on the currently known FRGs and IRGs and on their expression data from publicly available datasets. HCC patients with high CIFI values were found to have significantly shorter OSs, RFSs, DSSs, and PFIs than those with low CIFI values (Figures 3B,C, 4C-F). Moreover, the high-CIFI subgroup was found to have a significantly higher percentage of HCC patients with worse clinicopathological features, such as an advanced TNM stage, an advanced BCLC stage, cirrhosis, and vascular invasion ( Table 2). The patients in the high-CIFI subgroup showed a ferroptosis-suppressive status, as a set of well-defined SOFs was significantly upregulated in this subgroup (Figures 5A,B). A recent study showed that both SQSTM1 and NFE2L2 protect HCC cells from ferroptosis (Sun et al., 2016b); however, no consistent change was observed in these two genes between the highand low-CIFI subgroups, suggesting that the SQSTM1-Keap1-NFE2L2 pathway might not play a role in the CIFI-stratified subgroups. GPX4, a key regulator in suppressing ferroptosis (Hassannia et al., 2019;Seibt et al., 2019), was unexpectedly found to be down-regulated in the high-CIFI subgroup of the GSE14520 cohort, while its expression in the TCGA cohort was not different between the high and low-CIFI subgroups. This suggests that other imbalanced factors between the two subgroups, like hepatitis B virus (HBV) or hepatitis C virus (HCV) infection, or fatty liver disease, might play a predominant role in the expression of this gene. Indeed, other studies have found that GPX4 is over-expressed in HCV-related HCC patients and could be induced by HCV to increase virion infectivity (Guerriero et al., 2015;Brault et al., 2016). However, the HCV infection status of the HCC patients in these two cohorts was not available, and further studies may be required to address this question. In addition, two typical DOFs, namely NOX1 and ALOX15B, were found to be significantly downregulated in HCC patients with high CIFI values, while many other DOFs were not differentially expressed between the two subgroups (Supplementary Figure 1B). Interestingly, ACSL4 and ALOX5 were found to be significantly up-regulated in HCC patients with high CIFI values, even though they facilitate the execution of ferroptosis (Liu et al., 2015;Yuan et al., 2016;Doll et al., 2017). However, it is worth noting that both ACSL4 and ALOX5 are over-expressed in HCC as compared to normal livers, and promote the progression of the disease (Xu et al., 2011;Chen et al., 2020;Ndiaye et al., 2020;Wang et al., 2020). Thus, the over-expression of ACSL4 and ALOX5 contributes to the proliferation and progression of HCC, but renders the patients more susceptible to ferroptosis inducers. Indeed, researchers have found that some genes participate in multiple biological activities, some of which even seem contradictory. For instance, the well-known tumor suppressor TP53 plays an oncogenic role in HCC by inducing the P53 upregulated modulator of apoptosis (PUMA) (Kim et al., 2019). Nicotinamide phosphoribosyltransferase (NAMPT) drives tumor immune evasion by inducing PD-L1 expression in multiple types of tumors, but enhances the efficacy of immune checkpoint inhibitors (Lv et al., 2020). Further, GSEA revealed that genes in the low-CIFI subgroup were significantly enriched in biological oxidation, fatty acid metabolism, peroxisome, and glycine degradation ( Figure 5C). These biological processes are all critical for the execution of ferroptosis (Hassannia et al., 2019), suggesting a high sensitivity to ferroptosis in this subgroup. In addition, recent studies have found that ferroptosis inducers, such as erastin and sorafenib, could alleviate liver fibrosis (Sui et al., 2018;Wang L. et al., 2019;; thus, the ferroptosis-suppressive status observed in the high-CIFI subgroup might facilitate the progression of liver fibrosis and the development of cirrhosis. After all, the high-CIFI subgroup had a significantly higher percentage of HCC patients with cirrhosis ( Table 2).
Hepatocellular carcinoma patients with high CIFI values also exhibited immune-suppressive features. ssGSEA revealed that the CIFI value was positively correlated with the infiltration of MDSCs, which are able to promote immune escape and impair antitumor T-cell response in HCC (Chiu et al., 2017;Lin et al., 2017). On the contrary, eosinophil, which possesses anti-tumor activity toward tumors, including HCC (Kataoka et al., 2004;Hollande et al., 2019), was found to have a significantly negative association with the CIFI value in HCC ( Figure 6A). Unexpectedly, a positive association was also observed between CIFI value and CD4 T cells, which are regarded as cytotoxic to tumors (Jin et al., 2018). However, a study has shown that HBV-specific CD4 T cells are less cytotoxic to HCC cells and suppress the anti-tumor function of CD8 T cells (Meng et al., 2017). Indeed, 95.9% (211/220) of patients in the GSE14520 cohort were either active viral replication chronic carriers (AVR-CCs) or HBV chronic carriers (HBV-CCs) ( Table 1). Consequently, the positive relationship between CIFI and CD4 T cells might suggest an immunesuppression status instead of an immune-activation status. CIFI was also found to have a positive relationship with the fraction of CAFs in HCC ( Figure 6B). As one of the most abundant and critical components of the TME, CAFs contribute to immune evasion and immunotherapy failure, and promote the proliferation and invasion of tumors, including HCC, by secreting various growth factors and cytokines (Kubo et al., 2016;Chakravarthy et al., 2018;Liu T. et al., 2019). CIBERSORT analysis further indicated that HCC patients with high CIFI values had a significantly higher fraction of Tregs, which have been shown to be enriched in HCC and to inhibit IFN-gamma secretion and the cytotoxicity of CD8 + T cells (Yang et al., 2012;Langhans et al., 2019). Although less lytic than activated NK cells, resting NK cells are still cytotoxic and target cells like tumor cells, and could be converted to activated NK cells with certain stimuli (Bryceson et al., 2006;Lugini et al., 2012); however, the CIBERSORT analysis also revealed that the high-CIFI subgroup exhibited a lower infiltration of resting NK cells, while no difference was observed in the activated NK cells between the two subgroups ( Figure 6D). In addition, although the TMB was not found to be different between the two subgroups based on the calculation of gene mutation data, MATH analysis suggested that the high-CIFI subgroup had higher levels of tumor heterogeneity, which is generally a predictor of worse prognosis and is correlated with less immune response in various types of cancer (Mroz et al., 2015;Ma et al., 2019;McDonald et al., 2019). Considered together, these data suggest that high CIFI values might be correlated with immunosuppression in HCC.
In particular, it was observed that patients in the high-CIFI subgroup had a significantly higher frequency of TP53 mutation (41 vs. 24%). TP53 is a typical tumor suppressor, and its mutation leads to tumorigenesis and the progression of many types of tumors, including HCC (Leroy et al., 2014). However, TP53 mutation could either induce or suppress ferroptosis depending on the mutation site of the gene (Jiang et al., 2015;Jennis et al., 2016;Ou et al., 2016), and further studies are required to illustrate its exact impact on ferroptosis in HCC. Moreover, TP53 can activate an anti-tumor immune response via multiple mechanisms, like the down-regulation of immune-evading signals such as PD-L1, or the up-regulation of the NK cell ligands ULBP1 and ULBP2 (Textor et al., 2011;Cortez et al., 2016;Munoz-Fontela et al., 2016). HCC patients with TP53 WT have a significantly stronger local immuno-phenotype than those with mutant TP53 (Long et al., 2019). Consequently, the higher frequency of TP53 mutation observed in the HCC patients in the high-CIFI subgroup might contribute to the ferroptosis-and immune-suppressive phenotype of the subgroup.
A recent study predicted the prognosis of HCC patients by constructing a FRG signature from 60 FRGs . However, the results were not repeated and the signature was not compared with the model developed in the present work because the CARS gene could not be found in TCGA, GSE14520, or the Gene Cards website. 8 The authors might have mistaken CARS1 or CARS2 for CARS. In addition, the AUCs of the 10gene signature in the TCGA and ICGC cohorts were respectively 0.668 and 0.718 at 3 years; in contrast, those determined by the model in the present work were 0.72 and 0.7 at 3 years in the GSE14520 and TCGA cohorts, suggesting a slightly better predictability of the proposed model. In addition, the proposed signature was also found to have good predictability in longterm follow-up, as it achieved AUCs of 0.77 and 0.73 at 5 years in the GSE14520 and TCGA cohorts, respectively. Moreover, not only IRGs, but also FRGs (n = 283), were included in this work, and both suppressors and DOFs were considered for the construction of the predictive model. Therefore, the model might more comprehensively reflect the ferroptosis and immune status of HCC patients. 8 https://www.genecards.org/ The predictability of the proposed signature failed to be validated in the ICGC cohort because only a few dozen of the samples included the gene expression data of the selected genes. On the other hand, it should be noted that this work was limited because it was a retrospective study; thus, a further well-designed prospective analysis is necessary to validate the value of the developed model.
Because immunotherapy and targeted therapies, like sorafenib, function through ferroptosis, the novel comprehensive ferroptosis-immune status classifier developed in the present study suggests that personalized treatment should be applied in different subgroups of HCC patients. The high-CIFI subgroup represents a ferroptosis-and immune-suppressive phenotype, and might not benefit greatly from immunotherapy, targeted therapy, or a combined therapy like atezolizumab plus bevacizumab. On the other hand, the distinct features of this subgroup also imply that a ferroptosis inhibitor might synergize with immunotherapy and targeted therapy for a better treatment response.

CONCLUSION
In conclusion, a novel prognostic classifier based on ferroptosis and immune expression profiles in HCC was developed and validated. This classifier could be used for prognostic prediction and the selection of patients for immunotherapies and targeted therapies.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
ZS conceived and designed the study and drafted the manuscript. YL, XZ, and ZS carried out the data analysis. ZS and JL conducted the bioinformatics analysis with R software. YL and JT conducted the searching and reading of the existing literature. All authors revised the article and approved the final version to be published.