A Novel Ferroptosis-related Gene Signature for Overall Survival Prediction in Patients with Hepatocellular Carcinoma

1Liver cancer ranks sixth in terms of incidence among malignancies and is the fourth leading cause of tumor-related death worldwide. Hepatocellular carcinoma (HCC), the most prevalent type of primary liver cancer, is correlated with several well-known etiologies, including chronic HBV or HCV infection, alcohol abuse, nonalcoholic fatty liver disease, and exposure to dietary toxins such as aflatoxins. HCC is a highly heterogeneous disease that has been documented at interpatient, intertumoral and intratumoral level -. The overall survival of patients with HCC varies significantly across the world, with a 5-year survival rate of only 18% in the United States. The complex etiologic factors, along with the high-level heterogeneity of HCC, make prognostic prediction challenging. Furthermore, considering the limited treatment strategies for HCC, there is an additional need for the development of novel prognostic models.

2Ferroptosis is an iron-dependent form of regulated cell death (RCD) that is driven by the lethal accumulation of lipid peroxidation,. In recent years, the induction of ferroptosis has emerged as a promising therapeutic alternative to trigger cancer cell death, especially for malignancies that are resistant to traditional treatments,. Apart from ferroptosis-inducing agents, numerous genes have also been identified as modulators or markers of ferroptosis. Previous studies have reported that ferroptosis plays a vital role in HCC, and some genes, such as CISD1 and the polymorphism of the TP53 gene (S47 variant), are known to regulate ferroptosis negatively. On the other hand, other ferroptosis-related genes, such as Rb, NRF2 and MT1G, might protect HCC from sorafenib-induced ferroptosis. However, whether these ferroptosis-related genes are correlated with HCC patient prognosis remains largely unknown.

3In the present study, we first downloaded mRNA expression profiles and corresponding clinical data of HCC patients from public databases. Then, we constructed a prognostic multigene signature with ferroptosis-related differentially expressed genes (DEGs) in the TCGA cohort and validated it in the ICGC cohort. Finally, we further performed functional enrichment analysis to explore the underlying mechanisms.

4The level 3 RNA sequencing (RNA-seq) data and corresponding clinical information of 371 HCC patients were downloaded from the TCGA website up to November 15, 2019 (https://portal.gdc.cancer.gov/repository). The gene expression profiles were normalized using the scale method provided in the "limma" R package. RNA-seq data and clinical information of another 231 tumor samples were obtained from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). These samples were primarily derived from a Japanese population with HBV or HCV infection. Normalized read count values were used. The data from TCGA and ICGC are both publicly available. Thus, the present study was exempted from the approval of local ethics committees. The current research follows the TCGA and ICGC data access policies and publication guidelines.

5Then, 60 ferroptosis-related genes were retrieved from the previous literature,,, and are provided in Supplementary Table S1.

6The "limma" R package was used to identify the differentially expressed genes (DEGs) between tumor tissues and adjacent nontumorous tissues with a false discovery rate (FDR) < 0.05 in the TCGA cohort. Univariate Cox analysis of overall survival (OS) was performed to screen ferroptosis-related genes with prognostic values. P values were adjusted by Benjamini & Hochberg (BH) correction. An interaction network for the overlapping prognostic DEGs was generated by the STRING database (version 11.0). To minimize the risk of overfitting, the LASSO-penalized Cox regression analysis was applied to construct a prognostic model,. The LASSO algorithm was used for variable selection and shrinkage with the "glmnet" R package. The independent variable in the regression was the normalized expression matrix of candidate prognostic DEGs, and the response variables were overall survival and status of patients in the TCGA cohort. Penalty parameter (λ) for the model was determined by tenfold cross-validation following the minimum criteria (i.e. the value of λ corresponding to the lowest partial likelihood deviance). The risk scores of the patients were calculated according to the normalized expression level of each gene and its corresponding regression coefficients. The formula was established as follows: score= esum (each gene's expression × corresponding coefficient). The patients were stratified into high-risk and low-risk groups based on the median value of the risk score. Based on the expression of genes in the signature, PCA was carried out with the "prcomp" function of the "stats" R package. Besides, t-SNE were performed to explore the distribution of different groups using the "Rtsne" R package. For the survival analysis of each gene, the optimal cut-off expression value was determined by the "surv_cutpoint" function of the "survminer" R package. The "survivalROC" R package was used to conduct time‐dependent ROC curve analyses to evaluate the predictive power of the gene signature.

7The "clusterProfiler" R package was utilized to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the DEGs (|log2FC| ≥ 1, FDR < 0.05) between the high-risk and low-risk groups. P values were adjusted with the BH method. The infiltrating score of 16 immune cells and the activity of 13 immune-related pathways were calculated with single-sample gene set enrichment analysis (ssGSEA) in the "gsva" R package. The annotated gene set file is provided in Supplementary Table S2.

8Student's t-test was used to compare gene expression between tumor tissues and adjacent nontumorous tissues. Differences in proportions were compared by the Chi-squared test. Mann-Whitney test with P values adjusted by the BH method was used to compare the ssGSEA scores of immune cells or pathways between the high risk and low risk group. The OS between different groups was compared by Kaplan-Meier analysis with the log-rank test. Univariate and multivariate Cox regression analyses were implemented to identify independent predictors of OS. All statistical analyses were performed with R software (Version 3.5.3) or SPSS (Version 23.0). If not specified above, a P value less than 0.05 was considered statistically significant, and all P values were two-tailed.

9The flow chart of this study is shown in Fig. 1. A total of 365 HCC patients from the TCGA-LIHC cohort and 231 HCC patients from the ICGC (LIRI-JP) cohort were finally enrolled. The detailed clinical characteristics of these patients are summarized in Table 1.

10Most of the ferroptosis-related genes (49/60, 81.7%) were differentially expressed between tumor tissues and adjacent nontumorous tissues, and 27 of them were correlated with OS in the univariate Cox regression analysis (Fig. 2a). Unreasonably, HMOX1 was upregulated in tumor samples, but its expression predicted an excellent prognosis in the univariate Cox analysis, so it was excluded from further study. A total of 26 prognostic ferroptosis-related DEGs were preserved (all FDR <0.05, Fig.2b-c). The interaction network among these genes indicated that GPX4, G6PD and NQO1 were the hub genes (Fig.2d). The correlation between these genes is presented in Fig. 2e.

11LASSO Cox regression analysis was applied to establish a prognostic model using the expression profile of the 26 genes mentioned above. A 10-gene signature was identified based on the optimal value of λ (Fig. S1). Survival analyses, according to the optimal cut-off expression value of each gene, indicated that high expression of these genes all correlated with a poor prognosis (all adjusted P<0.05, Fig. S2). The risk score was calculated as follows: e (0.105 * expression level of SLC7A11 + 0.116 * expression level of G6PD + 0.106 * expression level of CISD1 + 0.076 * expression level of CARS + 0.077 * expression level of SLC1A5 + 0.092 * expression level of ACACA + 0.005 * expression level of ACSL3 + 0.006 * expression level of NQO1 + 0.087 * expression level of NFS1 + 0.135 * expression level of GPX4). The patients were stratified into a high-risk group (n=182) or a low-risk group (n=183) according to the median cut-off value (Fig. 3a). The higher risk group was found to be significantly associated with higher tumor grade, higher AFP, a higher rate of vascular invasion and advanced TNM stage in the TCGA cohort (Table 2). PCA and t-SNE analysis indicated the patients in different risk groups were distributed in two directions (Fig. 3b-c). As shown in Fig. 3d, patients with high risk had a higher probability of death earlier than those with low risk. Consistently, the Kaplan-Meier curve indicated that patients in the high-risk group had a significantly worse OS than their low-risk counterparts (Fig. 3e, P< 0.001). The predictive performance of the risk score for OS was evaluated by time-dependent ROC curves, and the area under the curve (AUC) reached 0.800 at 1 year, 0.690 at 2 years, and 0.668 at 3 years (Fig. 3f).

12Survival analyses of ten genes in the signature confirmed that these genes correlated with poor OS in the ICGC cohort except for CARS (all adjusted P< 0.05, Fig. S3). To test the robustness of the model constructed from the TCGA cohort, the patients from the ICGC cohort were also categorized into high- or low-risk groups by the median value calculated with the same formula as that from the TCGA cohort. The high risk group was also correlated with advanced TNM stage in the ICGC cohort (Table 2). Similar to the results obtained from the TCGA cohort, PCA and t-SNE analysis confirmed that patients in two subgroups were distributed in discrete directions (Fig. 4b-c). Likewise, patients in the high-risk group were more likely to encounter death earlier (Fig. 4d) and had a reduced survival time compared with those in the low-risk group (Fig. 4e, P= 0.001). Besides, the AUC of the 10-gene signature was 0.680 at 1 year, 0.690 at 2 years, and 0.718 at 3 years (Fig. 4f).

13Univariate and multivariate Cox regression analyses were carried out among the available variables to determine whether the risk score was an independent prognostic predictor for OS. In univariate Cox regression analyses, the risk score was significantly associated with OS in both the TCGA and the ICGC cohort (HR= 2.018, 95% CI = 1.409-2.891, P< 0.001; HR= 2.913, 95% CI = 1.489-5.699, P= 0.006, respectively) (Fig. 5a, b). After correction for other confounding factors, the risk score still proved to be an independent predictor for OS in the multivariate Cox regression analysis (TCGA cohort: HR =1.953, 95% CI = 1.356-2.812, P<0.001; ICGC cohort: HR=2.695, 95% CI = 1.357-5.351, P = 0.005; Fig. 5a, b).

14To elucidate the biological functions and pathways that were associated with the risk score, the DEGs between the high-risk and low-risk groups were used to perform GO enrichment and KEGG pathway analyses. As expected, DEGs were enriched in several iron-related molecular functions, such as ion channel activity and ion gated channel activity in both the TCGA and ICGC cohort (P. adjust < 0.05, Fig. 6a, c). Interestingly, the DEGs from the TCGA cohort were also obviously enriched in many immune-related biological processes (P. adjust < 0.05, Fig. 6a). Six immune-related biological processes or molecular functions were validated by the ICGC cohort, including chemokine-mediated signaling pathway, positive regulation of cytokine secretion, regulation of granulocyte macrophage colony-stimulating factor production, granulocyte macrophage colony-stimulating factor production, cytokine activity, and cytokine receptor binding (P. adjust < 0.05, Fig. 6c). KEGG pathway analyses also indicated that the cytokine-cytokine receptor interaction pathway was enriched in both cohorts (P. adjust < 0.05, Fig. 6b, d).

15To further explore the correlation between the risk score and immune status, we quantified the enrichment scores of diverse immune cell subpopulations, related functions or pathways with ssGSEA. Interestingly, contents of the antigen presentation process, including the score of aDCs, iDCs, APC co-stimulation, HLA and MHC class I, were significantly different between the low risk and high risk group in the TCGA cohort (all adjusted P< 0.05, Fig. 7a-b). The cytokine-cytokine receptor interaction that was enriched in the KEGG analyses had a higher score in the high risk group of the TCGA cohort (adjusted P< 0.05, Fig. 7b). Moreover, the score of type II IFN response, type I IFN response, and NK cells were lower in the high risk group, while the activity of checkpoint molecules, the scores of macrophages or Treg cells were just the opposite (adjusted P< 0.05, Fig. 7a-b). Comparisons in the ICGC cohort confirmed the differences of HLA, MHC class I, type II IFN response, checkpoint molecules, macrophages and Treg cells between two risk groups (adjusted P< 0.05, Fig. 7c-d). In particular, the scores of macrophages were the most statistically different between the two risk groups in both the TCGA and the ICGC cohort, which was consistent with the findings in the GO analysis.

16In the current study, we systematically investigated the expression of 60 ferroptosis-related genes in HCC tumor tissues and their associations with OS. A novel prognostic model integrating 10 ferroptosis-related genes was firstly constructed and validated in an external cohort. Functional analyses revealed that immune-related pathways were enriched.

17Although a few previous studies - have indicated that several genes might regulate drug-induced ferroptosis in HCC, their correlation with HCC patients' OS remains largely unknown. To our surprise, most of the ferroptosis-related genes (81.7%) were differentially expressed between tumor and adjacent nontumorous tissues, and more than half of them were correlated with OS in the univariate Cox regression analysis. These results significantly indicated the potential role of ferroptosis in HCC and the possibility of building a prognostic model with these ferroptosis-related genes.

18The prognostic model proposed in the present study was composed of 10 ferroptosis-related genes (ACACA, ACSL3, CISD1, CARS, G6PD, GPX4, NQO1, NFS1, SLC7A11, SLC1A5). These genes could be roughly classified into four categories, including iron metabolism (NFS1, CISD1), lipid metabolism (ACACA, ACSL3, GPX4), (anti)oxidant metabolism (CARS, NQO1, SLC7A11) and energy metabolism (G6PD, SLC1A5). NFS1, an enzyme involved in synthesizing iron-sulfur clusters using sulfur from cysteine, protects cells from ferroptosis in lung cancer. The genetic inhibition of CISD1 results in iron accumulation and subsequent oxidative injury in the mitochondria and thus contributes to erastin-induced ferroptosis in HCC cells. In terms of lipid metabolism, ACACA impacts the rate-limiting step in fatty acid synthesis and is a key regulator of tumor cell survival. Knockout of ACACA suppresses pharmacological agent-induced ferroptosis. ACSL3, which converts exogenous monounsaturated fatty acids (MUFAs) into fatty acyl-CoAs, is required for MUFA activation and promotes a ferroptosis-resistant cell status. GPX4 has been considered to be the primary enzyme that prevents ferroptosis for a long time due to its role in converting lipid hydroperoxides into nontoxic lipid alcohols. For (anti)oxidant metabolism, knockdown of CARS inhibits erastin-induced ferroptosis by preventing the induction of lipid reactive oxygen in fibrosarcoma cells. In contrast, knockdown of NQO1 enhances erastin and sorafenib-induced ferroptosis in HCC cells. Similarly, knockdown of SLC7A11, a subunit of system Xc to import cystine in the cell, sensitized fibrosarcoma cells to erastin-induced death. G6PD and SLC1A5 are two ferroptosis regulators involved in energy metabolism. G6PD, which is involved in the pentose phosphate pathway, has been reported to prevent erastin-induced ferroptosis when it was knocked down in non-small cell lung cancer cells. Inhibition or knockdown of SLC1A5 also suppresses ferroptosis. In summary, six of the genes (NFS1, CISD1, ACSL3, NQO1, SLC7A11, GPX4) in the prognostic model have been reported to protect cells from ferroptosis, while the remaining four genes (ACACA, CARS, G6PD, SLC1A5) are the opposite. However, these genes were all upregulated in HCC tumor tissue and were associated with poor prognosis in the current study. Whether these genes play a role in HCC patients' prognosis by influencing the process of ferroptosis remains to be elucidated, since few related studies on these genes except for CISD1 and NQO1 have been reported.

19Although the mechanisms underlying tumor susceptibility to ferroptosis have been an intense area of research in the past few years, the potential modulation between tumor immunity and ferroptosis remains elusive. Based on the DEGs between different risk groups, we performed GO analyses and unexpectedly discovered that many immune-related biological processes and pathways were enriched. It is reasonable to assume that ferroptosis may have a close connection with tumor immunity. Interestingly, the contents of the antigen presentation process were significantly different between the low risk and high risk group in this study. One possible speculation is that ferroptotic cells release distinct signals, such as lipid mediators, to attract antigen-presenting cells (APCs) to the site of ferroptotically dying cells. Besides, the high risk groups in both the TCGA and the ICGC cohort have higher fractions of macrophages and Treg cells. Previous studies have demonstrated that increased tumor-associated macrophages [30, 31]or Treg cells [31, 32]are related to poor prognosis in HCC patients due to their role in immune invasion. Moreover, higher risk scores correlated with impaired antitumor immunity, including the activity of the type II IFN response and type I IFN response as well as the fractions of NK cells. Therefore, attenuated antitumor immunity in patients with high risk may be an explanation for their poor prognosis.

20There are several limitations of this study. First, our prognostic model was both constructed and validated with retrospective data from public databases. More prospective real-world data are warranted to verify its clinical utility. Second, the intrinsic weakness of merely considering a single hallmark to build a prognostic model was inevitable, as many prominent prognostic genes in HCC might have been excluded. In addition, it should be emphasized that the links between the risk score and immune activity have not yet been experimentally addressed.

21In summary, our study defined a novel prognostic model of 10 ferroptosis-related genes. This model proved to be independently associated with OS in both the derivation and validation cohorts, providing insight into the prediction of HCC prognosis. The underlying mechanisms between ferroptosis-related genes and tumor immunity in HCC remain poorly understood and warrant further investigation.