Integrated Analysis of Immunity- and Ferroptosis-Related Biomarker Signatures to Improve the Prognosis Prediction of Hepatocellular Carcinoma

Background Hepatocellular carcinoma (HCC) is a common malignant tumor with high mortality and poor prognoses around the world. Ferroptosis is a new form of cell death, and some studies have found that it is related to cancer immunotherapy. The aim of our research was to find immunity- and ferroptosis-related biomarkers to improve the treatment and prognosis of HCC by bioinformatics analysis. Methods First, we obtained the original RNA sequencing (RNA-seq) expression data and corresponding clinical data of HCC from The Cancer Genome Atlas (TGCA) database and performed differential analysis. Second, we used immunity- and ferroptosis-related differentially expressed genes (DEGs) to perform a computational difference algorithm and Cox regression analysis. Third, we explored the potential molecular mechanisms and properties of immunity- and ferroptosis-related DEGs by computational biology and performed a new prognostic index based on immunity- and ferroptosis-related DEGs by multivariable Cox analysis. Finally, we used HCC data from International Cancer Genome Consortium (ICGC) data to perform validation. Results We obtained 31 immunity (p < 0.001)- and 14 ferroptosis (p < 0.05)-related DEGs correlated with overall survival (OS) in the univariate Cox regression analysis. Then, we screened five immunity- and two ferroptosis-related DEGs (HSPA4, ISG20L2, NRAS, IL17D, NDRG1, ACSL4, and G6PD) to establish a predictive model by multivariate Cox regression analysis. Receiver operating characteristic (ROC) and Kaplan–Meier (K–M) analyses demonstrated a good performance of the seven-biomarker signature. Functional enrichment analysis including Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed that the seven-biomarker signature was mainly associated with HCC-related biological processes such as nuclear division and the cell cycle, and the immune status was different between the two risk groups. Conclusion Our results suggest that this specific seven-biomarker signature may be clinically useful in the prediction of HCC prognoses beyond conventional clinicopathological factors. Moreover, it also brings us new insights into the molecular mechanisms of HCC.


INTRODUCTION
Liver cancer is a common malignant tumor around the world. Liver cancer is the seventh most common cancer globally according to global cancer data in 2017, with 953,000 new cases diagnosed and 819,000 deaths (Fitzmaurice et al., 2019). Hepatocellular carcinoma (HCC) is the most common because HCC accounts for approximately 90% of primary liver cancer (Llovet et al., 2016). Risk factors for HCC mainly include cirrhosis (chronic liver damage caused by inflammation and fibrosis), hepatitis B virus (HBV) infection, hepatitis C virus (HCV) infection, alcohol abuse, and non-alcoholic fatty liver disease (NAFLD) (Fujiwara et al., 2018). HCC is very malignant, and 70% of patients undergoing surgery experience tumor recurrence within 5 years (Nakagawa et al., 2016). In addition, HCC is highly complex and heterogeneous, so many molecular targeted anticancer agents are not effective for some patients, with some even showing resistance (Cancer Genome Atlas Research Network, 2017). The prognosis of HCC is very poor, with a 3 years survival rate of 12.7% and a 5 years survival rate of 20% (Giannini et al., 2015;Goossens et al., 2015;Yu et al., 2017). Currently, the diagnosis and treatment of HCC are not satisfactory. Specific biomarkers play an important role in the early screening, diagnosis, treatment option selection, and prognosis of HCC. In our study, we explored some immunityand ferroptosis-related biomarkers to understand how they affect the pathogenesis and prognosis of HCC. We hope that these biomarkers will be helpful for the diagnosis, treatment, and prognosis of HCC.
Ferroptosis is an iron-dependent form of non-apoptotic cell death in the presence of oxidized polyunsaturated fatty acids (PUFAs), redox-active iron, and compromised lipid peroxide repair (Perez et al., 2020). Ferroptosis has been demonstrated in many diseases, such as cancer (Liang et al., 2019;Perez et al., 2020). It also plays a very vital role in digestive-system neoplasms such as gastric cancer, pancreatic cancer, colorectal cancer, and, especially, HCC (Nie et al., 2018). Research shows that the p62-Keap1-NRF2 pathway is related to the ferroptosis of liver cancer cells, and the retinoblastoma (RB) protein plays a role in liver tumorigenesis and ferroptosis (Viatour et al., 2011;Sun et al., 2016;Nie et al., 2018). Therefore, it is important to research ferroptosis-related biomarkers in the prognosis and treatment of HCC patients.
In total, 80% of HCC patients in an advanced stage lost opportunities for surgery and ablation, but systemic treatments for HCC are limited (Liu and Qin, 2019). In recent years, immunotherapy has been given more attention, and clinical trials and animal experiments have proven that immunotherapy plays a role in the treatment of HCC patients (Zongyi and Xiaowu, 2020). Immunotherapy approaches, including vaccines, immune checkpoint blockade, and adoptive cell transfer (ACT), have been proven safe and effective for HCC treatment (Li et al., 2015). Therefore, it is important to explore immunity-related biomarkers for immunotherapy treatment of HCC.
In summary, we find that ferroptosis in HCC must be related to immunity. Lang et al. (2019) found that radiotherapy and immunotherapy can promote ferroptosis via SLC7A11. Ubellacker et al. (2020) proved that metastasizing melanoma cells from the lymph nodes are resistant to ferroptosis. Jiang et al. (2020) designed sulfasalazine (SAS)−loaded mesoporous magnetic nanoparticles (Fe 3 O 4 ) and platelet (PLT) membrane camouflage (Fe 3 O 4 -SAS @ PLT), which can mediate ferroptosis with immunotherapy and can be expected to provide great potential in the clinical treatment of tumor metastasis. Ruiz-de-Angulo et al. (2020) developed iron oxide-loaded nanovaccines (IONVs) that can enhance the combination immunotherapy and immunotherapy-promoted tumor ferroptosis. We find that exploring the clinical relevance and prognostic significance of immunity-and ferroptosis-related biomarkers is helpful for ferroptosis immunotherapy in HCC.
In this study, we used TCGA database to analyze the mRNA expression profiles to find immunity-and ferroptosis-related differentially expressed genes (DEGs) for the prognosis of HCC. Furthermore, functional analysis of potential immunityand ferroptosis-related DEGs is helpful for understanding their roles in the occurrence and development of HCC. Finally, we also validated our results in the International Cancer Genome Consortium (ICGC) cohort. Therefore, this study provides a good prognostic risk model for HCC patients and some insights into the occurrence and development of HCC.

Patients and Datasets
The RNA sequencing (RNA-seq) data and corresponding clinical information of 374 HCC samples and 50 normal liver samples were downloaded from TCGA database 1 on September 10, 2020. In addition, The RNA-seq data and corresponding clinical information of 231 liver cancer samples were obtained from the ICGC database 2 on September 10, 2020. From the Immunology Database and Analysis Portal (ImmPort), a list of 2,483 immunity-related genes was obtained. ImmPort is a database that updates immunology data accurately and in a timely manner. Data from ImmPort are a powerful foundation of immunology research. In addition, it provides a list of immunityrelated genes for use in cancer research. These genes were proven to participate in the process of immune activity. Additionally, 60 ferroptosis-related genes were obtained by summarizing previous literature.

Differential Gene Analysis
By comparing HCC tissues to normal tissues and using the limma package from Bioconductor in R software (version 4.0.2), DEGs were identified with the criteria for screening the DEGs of a false discovery rate (FDR) < 0.05 and a | log2FoldChange| > 1. Then, immunity-and ferroptosis-related DEGs were extracted from the DEGs.     Frontiers in Genetics | www.frontiersin.org

Construction and Validation of the Prognostic Immunity-and Ferroptosis-Related Differentially Expressed Gene Signature
The association between immunity-and ferroptosis-related DEGs and patient survival was evaluated by univariate Cox regression analysis using the survival R package in R. Immunityrelated DEGs with a p < 0.001 and ferroptosis-related DEGs with a p < 0.05 were considered candidate variables in a univariate Cox regression analysis and entered into a stepwise multivariate Cox regression analysis tested by Akaike Information Criterion (AIC, assessing the goodness of fit of a statistical model) to identify covariates with independent prognostic values for patient survival. Based on the median risk score, HCC patients were divided into high-and low-risk groups. The Kaplan-Meier (K-M) survival curves for cases predicted with low and high risk were generated, respectively. Then, ROC curve analysis was performed to test the sensitivity and specificity of the prognostic risk score model for overall survival (OS) in HCC patients. The area under the ROC curve (AUC) was derived as reported previously. Then, we performed an independent prognostic analysis to evaluate whether clinical parameters are independent prognostic factors for OS. Based on prognostic immunity-and ferroptosis-related DEGs, we performed principal component analysis (PCA) by the "prcomp" function of the "stats" R package and used t-distributed stochastic neighbor embedding (t-SNE) with the "Rtsne" R package. To verify these results, we used ICGC data.

Construction of the Immunity-and Ferroptosis-Related Differentially Expressed Gene Regulatory and Protein-Protein Interaction Network
We explored prognostic immunity-and ferroptosis-related DEG regulatory mechanisms, extracted prognostic immunityrelated DEGs to construct the regulatory network of immunity-and ferroptosis-related DEGs, and used Cytoscape software (version 3.7.2) to display the immunity-and ferroptosis-related DEG regulatory network. To explore the interactions between these genes, we constructed a protein-protein interaction (PPI) network based on data gleaned from the Retrieval of Interacting Genes (STRING) online database 3 .

Functional Enrichment Analysis
Based on the median risk score, HCC patients were divided into the high-and low-risk groups. We applied the "limma" R package to analyze the correlations of DEGs between the high-and low-risk groups with FDR < 0.05 and | log2FoldChange| > 1 in TCGA and ICGC cohorts, respectively. Then, we used riskrelated DEGs to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses by the "clusterProfiler" R package. Finally, for differential infiltrating score analysis between the highand low-risk groups, we determined the infiltrating scores of immune cells and immune-related functions for samples by single-sample gene set enrichment analysis (ssGSEA).

Identification of Immunity-and Ferroptosis-Related Genes
Using the limma package of R language for DEGs (FDR < 0.05 and | log2FoldChange | > 1) analysis, we obtained 7,768 DEGs in the 374 HCC samples and 50 normal liver samples from TCGA database. Then, we extracted 335 immunity-related DEGs from 2,483 entries of the ImmPort database and 26 ferroptosis-related DEGs from 60 ferroptosis-related genes in literature reporting ( Table 1). The flow diagram of this study is shown (Figure 1).

Identification of Prognostic Immunityand Ferroptosis-Related Differentially Expressed Genes
To explore the immunity-and ferroptosis-related DEG correlations with the OS of HCC patients, we obtained 31 immunity (p < 0.001)-and 14 ferroptosis (p < 0.05)-related DEGs by univariate Cox regression analysis (Tables 2, 3). According to the forest plot of hazard ratios (HRs), most of these DEGs were risk factors for poor prognoses in HCC patients. Then, we used these DEGs to perform multivariate Cox regression analysis. Finally, we identified five immunityand two ferroptosis-related DEGs (HSPA4, ISG20L2, NRAS, IL17D, NDRG1, ACSL4, and G6PD) to establish a predictive model ( Table 4). K-M survival curves outcomes based on median risk score values show that the predicted survival time of the low-risk group was significantly longer than that of the high-risk group, p < 0.001 (Figure 2A). We used the time-dependent ROC curve to assess the ability of the risk score to predict survival rates ( Figure 2B). The results showed that the AUC of the five immunity-and two ferroptosis-related DEGs in the prognostic model was 0.837 and demonstrated that the risk score model had a stable performance. HCC patients were classified into highand low-risk groups according to respective median risk scores (Figures 2C,D). By PCA and t-SNE analyses, we also indicated that HCC patients in different risk groups were distributed in two directions (Figures 2E,F).

Construction of the Immunity-and Ferroptosis-Related Differentially Expressed Gene Regulatory and Protein-Protein Interaction Network
First, we constructed a regulatory network about prognostic 31 immunity-and 14 ferroptosis-related DEGs ( Figure 3A). Then, using the STRING online platform, we established a PPI network based on these DEGs ( Figure 3B). Both networks indicated that immunity and ferroptosis have potential molecular mechanisms.

Validation of the Seven-Biomarker Signature Using the International Cancer Genome Consortium Data
To verify the reliability of the model from TCGA data, we selected seven biomarkers to perform multivariate Cox regression analysis using ICGC data. Likewise, compared to the median risk score values, the K-M survival curve outcomes show similar results ( Figure 4A). In addition, the time-dependent ROC curve showed that the AUCs of the seven-biomarker prognostic model was 0.787 ( Figure 4B). It also indicated that the risk score model was very robust. HCC patients with ICGC data were also categorized into high-and low-risk groups according to the respective median risk scores (Figures 4C,D). PCA and t-SNE analyses demonstrated that patients in the two groups were distributed in discrete directions (Figures 4E,F).

Independent Prognostic Analysis
We performed univariate and multivariate Cox regression analyses to evaluate whether clinical parameters (including gender, age, stage, and grade) and the risk score are independent prognostic factors of OS. Then, in both TCGA and ICGC data, we found that the risk score and stage were independent prognostic predictors for OS in the univariate and multivariate Cox regression analyses (Figures 5A-D).

Functional Enrichment Analyses
We performed GO and KEGG functional enrichment analyses on risk-related DEGs to investigate the potential functions of the seven prognostic biomarkers. The results indicated that the seven prognostic biomarkers mainly focused on nuclear division and mitotic nuclear division in both TCGA and ICGC data (Figures 6A,B). KEGG functional enrichment analysis suggested that the seven prognostic biomarkers were mainly related to the cell cycle, the metabolism of xenobiotics by cytochrome P450, extracellular matrix (ECM)-receptor interactions, etc. (Figures 6C,D). To further explore relationship between the HCC prognosis and immune status, we quantified the infiltrating scores of immune cell-and immunity-related functions with ssGSEA. The correlations between ssGSEA scores and different risk groups showed that the scores of iDCs, Macrophages, NKcells, Th2-cells, Treg, APC costimulation, Type I IFN Response, and Type II IFN Response were significantly different between the low-and high-risk groups in both TCGA and ICGC cohorts (Figures 7A-D).

DISCUSSION
Although the current treatment and diagnosis of HCC have improved compared with the past, the mortality and incidence of HCC are still high. Colorectal cancer, stomach cancer, and HCC have the highest mortality rates, which have surpassed those of lung cancer (Chi et al., 2019). The prognosis is poor. It is very important to predict the prognosis of HCC and give corresponding treatments in time. The OS rate of HCC is still very low, and the prognosis is poor. The quality of life with advanced HCC is also poor (Anwanwan et al., 2020). Therefore, it is important to have deeper insight into the pathological mechanisms of the occurrence and development of HCC. In addition, it is vital to find new biomarkers and targets that are more sensitive and effective for the early diagnosis, treatment, and prognosis of HCC.
At present, many researchers have proven that ferroptosis is related to immunity, and immunotherapy and ferroptosis are both current research hotspots. However, it is rare to construct prognostic models of immunity-and ferroptosisrelated genes and to predict their correlations in HCC. In our study, we systematically investigated the differential expression of immunity-and ferroptosis-related genes in HCC and performed a survival analysis. A novel prognostic model integrating five immunity-and two ferroptosis-related DEGs was first constructed and validated in an external cohort. Functional analyses show that immune cells and immunity-related functions and pathways were enriched.
We established the predictive model including five immunityand two ferroptosis-related DEGs (HSPA4, ISG20L2, NRAS, IL17D, NDRG1, ACSL4, and G6PD). At present, the full name of HSPA4/Hsp70 is Heat Shock Protein Family A Member 4, which is a protein-coding gene. Chen et al. (2009) showed that the FOLFOX4 (5-fluorouracil, leucovorin, and oxaliplatin) regimen has advantages in colorectal cancer patients who have unresectable liver metastasis with lower HSP70 expression over those with higher HSP70 expression. Zhu et al. (2019) proved that ISG20L2 was related to the HCC prognosis. Concerning NRAS, Dietrich et al. (2019) found that it is a prognostic biomarker and contributes to sorafenib resistance in HCC. O'Sullivan et al. (2014) proved that IL17D was a novel target for the immunotherapy of tumors. NDRG1 is a hypoxia-inducible protein, which is related to the progression of human cancers, and is induced by hypoxia in HCCs (Sibold et al., 2007). Long non-coding RNA CCAT2 promotes proliferation and metastasis through upregulation of NDRG1 in HCC . The FOXQ1/NDRG1 axis can drive HCC initiation; thus, NDRG1 was suggested as a new potential therapeutic target for HCC (Luo et al., 2018). ACSL4 was overexpressed and served as an independent adverse prognostic index in HCC (Sun and Xu, 2017). Lu et al. (2018) found that overexpression of G6PD was correlated with epithelial-mesenchymal transition, which contributes to the migration and invasion of HCC.
At present, cancer treatment has entered an age of immunity and iron (Buttner et al., 2016;Tarangelo and Dixon, 2016). For ferroptosis, nanoparticles modulate iron and ROS levels to induce ferroptosis, which could provide a new strategy for cancer treatment (Tarangelo and Dixon, 2016). Immunotherapy has become a new standard of treatment for advanced HCC around the world (Waidmann, 2018). CD8 + T cells release interferon (IFN)γ, and IFNγ can regulate the lipid peroxidation and ferroptosis pathways in tumors . These findings mean that ferroptosis can play a role in the immunotherapy of tumors.
There are several limitations to our study. First, the number of normal samples from TCGA database was relatively small. Second, the present study only included a database mining design, without validation using fresh samples and prospective experimental research. Therefore, we will collect fresh samples and further prove our conclusions through experiments.

CONCLUSION
In conclusion, we identified seven biomarkers associated with the prognosis of HCC patients. The prediction of the seven biomarkers' functions provided further insight into the roles of the seven biomarkers in the immunotherapy and ferroptosis therapy of HCC. We also proved that the seven-biomarker signature can well predict the survival times of HCC patients. This signature has many potential prognostic and therapeutic implications for HCC patient management.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: http://portal.gdc.cancer.gov/repository.

AUTHOR CONTRIBUTIONS
YZ conceived and designed the study. XD conducted the analysis and wrote the manuscript. Both authors read and approved the final manuscript.