Construction and validation of a prognostic model for hepatocellular carcinoma: Inflammatory ferroptosis and mitochondrial metabolism indicate a poor prognosis

Abstract

Background:

An increasing number of innovations have been discovered for treating hepatocellular carcinoma (HCC or commonly called HCC) therapy, Ferroptosis and mitochondrial metabolism are essential mechanisms of cell death. These pathways may act as functional molecular biomarkers that could have important clinical significance for determining individual differences and the prognosis of HCC. The aim of this study was to construct a stable and reliable comprehensive model of genetic features and clinical factors associated with HCC prognosis.

Methods:

In this study, we used RNA-sequencing (fragments per kilobase of exon model per million reads mapped value) data from the Cancer Genome Atlas (TCGA) database to establish a prognostic model. We enrolled 104 patients for further validation. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes enrichment analyses (KEGG) analysis were used for the functional study of differentially expressed genes. Pan-cancer analysis was performed to evaluate the function of the Differentially Expressed Genes (DEGs). Thirteen genes were identified by univariate and least absolute contraction and selection operation (LASSO) Cox regression analysis. The prognostic model was visualized using a nomogram.

Results:

We found that eight genes, namely EZH2, GRPEL2, PIGU, PPM1G, SF3B4, TUBG1, TXNRD1 and NDRG1, were hub genes for HCC and differentially expressed in most types of cancer. EZH2, GRPEL2 and NDRG1 may indicate a poor prognosis of HCC as verified by tissue samples. Furthermore, a gene set variation analysis algorithm was created to analyze the relationship between these eight genes and oxidative phosphorylation, mitophagy, and FeS-containing proteins, and it showed that ferroptosis might affect inflammatory-related pathways in HCC.

Conclusion:

EZH2, GRPEL2, NDRG1, and the clinical factor of tumor size, were included in a nomogram for visualizing a prognostic model of HCC. This nomogram based on a functional study and verification by clinical samples, shows a reliable performance of patients with HCC.

Introduction

Liver hepatocellular carcinoma (HCC) is the third leading cause of cancer-related death worldwide (1). Most patients with HCC have a history of chronic hepatitis B virus (HBV) infection. Although HBV has been effectively controlled in China, morbidity and mortality are currently rising (2). Nevertheless, HCC is the sixth most common cancer worldwide (3). Several innovative methods have been developed for the treatment of HCC. However, the effectiveness of these treatments markedly varies between individuals. The main reason for this variation is a lack of therapeutic evaluation of efficacy. Therefore, more effective biomarkers and/or prognostic models are requied. Effective biomarkers should be able to predict the prognosis of HCC and participate in specific functions of tumor biology. In particular, biomarkers described for specific anti-tumor phenotypes may be the most meaningful (4, 5).

Intracellular iron is a crucial target for the induction of cell death because it mediates ferroptosis and mitochondrial metabolism (6). Ferroptosis is a non-apoptotic regulated form of cell death that has gradually attracted attention in cancer pathogenesis research (7). Ferroptosis induces abnormal mitochondrial metabolism and reactive oxygen species (ROS) accumulation, and it may promote anti-tumor mechanisms in the microenvironment (8, 9). These finding suggests that ferroptosis is an important mechanism contributing to the development of the tumor microenvironment. Some researchers have suggested that ferroptosis or metabolic-related gene signatures could provide relevant information on the prognosis of HCC (10, 11). The unique tumor microenvironment in HCC includes considerable immune infiltration because of chronic inflammation (12). Therefore, a prognostic model for HCC in the context of inflammation in mainland China is required.

In this study, a bioinformatics analysis method was used to analyze a data set from the Cancer Genome Atlas (TCGA) database to determine the differentially expressed genes (DEGs) in HCC compared with healthy individuals. A prognostic model based on the hub genes was constructed using least absolute contraction and selection operation (LASSO) regression. To determine the validity and stability of this model, we used clinical tissue samples to verify the hub genes. To verify the function of these critical genes in tumors, we performed a series of functional analyses. Using these results, we hoped to identify the most stable and reliable gene signatures and clinical factors that contribute to the progression of HCC.

Materials and methods

Data sources

Our analysis results are based on the omics data set generated by the TCGA (http://cancergenome.nih.gov/). A total of 33 types of cancer were enrolled from TCGA. All data sets include Count data, Transcripts Per Kilobase Million (TPM) data, RNA-seq, Copy number variation (CNV), Mutation and Clinical Information, were downloaded from UCSC XENA (https://xenabrowser.net/). We enrolled a total of 50 normal liver tissues and 374 HCC in the TCGA database. We used the limma package for R to analyze the differential expression between HCC and normal liver tissues.

Clustering differentially expressed genes for HCC

To identify the Fold Changes in gene expression in each cancer data set, we used the Deseq2 package for R software to identify the DEGs using TCGA count data. We consider adjusting P value < 0.05 and Log|FC| > 1 as DEGs. Meanwhile, to determine the DEGs related to the prognosis of HCC, we used univariate COX regression, which is analyzed in R using survival package to screen the genes that affect the prognosis of HCC. HCC-specific genes are considered meaningful to both differential expression and prognostic analysis.

GO and KEGG analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes Enrichment Analyses were functional enrichment. Go analysis is a bioinformatics tool including three functional analyses, which are cellular components (CC), molecular function (MF), and biological pathways (BP) (13). KEGG stores a multigene signature pathway to which hub genes belong (14). We conducted a functional enrichment analysis to study the underlying mechanism of DEGs-based clustering at different levels. We first converted the gene names of the differential genes into ENTREZ IDs, and further, enriched using the enrichment function in the cluterprofiler. For the P value of the results of the analysis, the BH algorithm is used to correct it. Finally, we select the corrected P value of < 0.05 as the meaningful enrichment result

LASSO Cox regression analysis and visualization

We first performed COX regression on differentially expressed genes to screen genes related to the prognosis of HCC. The LASSO regression analysis model eliminates the collinearity problem between prognostic-related genes (15). Finally, we obtained the characteristic genes that affect HCC. We used a nomogram, which constructed specific models related to HCC for clinical application to show the clinical model. Lasso returns to the Glmnet package using the R language. We used the Rms package to construct the nomogram graph.

Patients and tissue samples collection

A total of 104 HCC patients were enrolled in The Cancer Hospital of the University of Chinese Academy of Sciences(Zhejiang Cancer Hospital) from 2011 to 2020. Each patient was collected one pair of normal and cancer tissue within 30 minutes after surgical resection. Normal tissue was defined as greater than 1cm beyond the tumor margin. A frozen section was identified by senior pathologists, who were more than 15 years of practitioners, to exclude tumor lesions under the microscope. HCC lesions, which senior pathologists diagnosed, were acquired about 1cm inside the tumor margin. All patients were pathologically confirmed as Hepatocellular carcinoma. The clinical stage of the tumor was determined according to the Cancer Staging Manual of the American Joint Committee on Cancer(version 8, 2017), while the Edmondson Steiner classification defined the tumor pathological differentiation stage. All 104 patients were enrolled, including 19 females and 85 males. There were 58 patients older than 60 years old and another 46, respectively. The follow-up ended in March 2021 or at death. HCV infections were excluded in our study.

Immunohistochemistry stain

All 104 pairs of lesions were paraffin-fixed tissue. All slice tissue was deparaffinized, rehydrated, and then heated in citric acid buffer (0.01M) at 105°C for 10 min. After antigen recovery, the sliced tissue was blocked with 3% hydrogen peroxide solution and bovine serum albumin. Subsequently, the slides were incubated with rabbit antibodies, including anti-EZH2, anti-GRPEL2, anti-PIGU, anti-PPM1G, anti-SF3B4, anti-TCOF1, anti-TUBG1, anti-TXNRD1, anti-MYCN, anti-NDRG1, anti-SQSTM1, anti-UCK2, at 4°C overnight, then incubated with horseradish peroxidase (HRP)-conjugated secondary antibody at room temperature for 20 min. Finally, the slides were counterstained with hematoxylin, dehydrated in graded alcohol and xylene, cleared, and mounted. The antibody or catalogue and dilutions are listed in Supplementary Table 1.

IHC quantification

All slices were measured by Pannoramic DESK II DW Digital Slide Scanner(Coherent Scientific Ltd.). After the tissue section was set on the scanner, the section gradually moved under the scanner lens. While moving, all tissue information would be scanned and recorded to form a folder containing all image details. The folder was opened by using CaseViewer2.2 software. It can be magnified at any multiple of 1x-400x for observation. Use the TMA plug-in in the Quant Center2.1 analysis software to set the diameter of the chip organization point and the number of rows and columns. The software will automatically generate the number. Use the PatternQuant module in the Quant Center2.1 analysis software to distinguish the brown area (including positive) from the blue part. The two parts of the Mask Area are used as the tissue area; add HistoQunant is merged in the next layer of the brown area. Finally, the area and grayscale of the positive part are calculated based on the optical density.

Pan-cancer analysis for DEGs

To determine the role of the above 13 DEGs, we used the genes to perform pan-cancer analysis to observe the expression, mutation frequency, and copy number change frequency of these genes in pan-cancer. All TCGA pan-carcina data can be downloaded from UCSC XENA. In gene mutation analysis, we calculated the proportion of mutations in each gene in different cancer samples to compare which gene was more likely to mutate. In copy number analysis, the proportion of each gene that increases or decreases in copy number is calculated. We observed changes in gene copy number using the increase in copy number and the decrease in ratio. Finally, we used Deseq2 to calculate the differential expression of each gene. The differential expression of core genes in multiple cancers and the effect of gene expression on the prognosis of different tumors were observed.We used TCGA+GTEx and GSE101685 to validate the core genes. TCGA+GTEx contains 369 cancers and 160 normals, while GSE101685 includes 24 cancer samples and eight normal samples.

Identification of hub genes and functional enrichment analysis

The Gene Set Variation Analysis(GSVA) algorithm is an algorithm that obtains phenotypic scores based on the amount of gene expression associated with different phenotypes. We first collected genes associated with each immune cell. Further, based on the expression of these genes, each immune cell score is obtained. First, we used the Gene Set Variation Analysis(GSVA) to calculate the impact of 24 types of immune cells on HCC. The relationship between core genes and immune cell scores was analyzed. Pearson correlation analysis was further used for the relationship between hub genes and immune cells in HCC. To investigate the relationship between risk score and immune checkpoints, we extracted the expression of 30 immune checkpoints, including the B7-CD28 family (CD274, CD276, CTLA4, HHLA2, ICOS, ICOSLG, PDCD1, PDCD1LG2, TMIGD2, VTCN1), TNF superfamily (BTLA, CD27, CD40LG, CD40, CD70, TNFRSF18, TNFRSF4, TNFRSF9, TNFSF14, TNFSF4, TNFSF9), and another immune checkpoint (HAVCR2, IDO1, LAG3, FGL1, ENTPD1, NT5E, SIGLEC15, VSIR, NCR3). Furthermore, we measured the coexpression relationship between DEGs and Ferroptosis. As all we know, cell death is affected by mitochondria dysfunction generally. Therefore, we downloaded the data from MITOCARTA 3.0 to evaluate mitochondria-related pathways. The GSEA algorithm was used to analyze the effect of DEGs on OXPHOS, Mitophagy and Fe.S containing proteins.

Identification of hub genes and drug sensitivity and regulatory mechanisms

We assessed the drug sensitivity of each patient using pRRophetic (R package, University of Minnesota System), and further used correlation analysis to relate the relationship between prognostic models and tumor drug sensitivity. pOSTAR3 is a database used to predict post-transcriptional regulation of genes, and to understand the transcriptional regulatory relationships of three iron death-associated genes, we constructed a post-transcriptional regulatory network using POSTAR3. GSE106988 is a dataset of RNA-seq performed after knockdown of NDRG1. Using this database. We analyzed the potential functions of NDRG1 and used cytoscape software for the visualization of protein pathway construction.

Statistical analysis

R (version 3.6.2, www.r-project.org) and relative packages for R were used for bioinformatics and statistical analysis. SPSS (Version 25.0) was used for statistical analysis. Visualization mainly used the R package ggplot2(version 3.3.3). The chi-square test and Fisher’s exact test were used to evaluate the statistical significance of the relationship between DEGs expression and clinicopathological parameters. The Student’s unpaired t-test is used to express customarily distributed DEGs; otherwise, the Mann-Whitney U test is used. The Wilcoxon signed-rank test was used to analyze the difference between HCC and adjacent non-HCC lesion distribution. Survival curves were assessed by the Kaplan-Meier method and Cox proportional hazards regression model, and a log-rank test analyzed the differences. Before performing the Kaplan-Meier analysis, the Best Separation method was used to cut gene expression into high and low groups. P < 0.05 was considered as statistical significance. After building the model, we evaluated our model by C-index and decision curve.

Results

Identification of DEGs from the TCGA

The entire analysis flow of the manuscript results is shown in Figure 1. DESeq2 was performed for DEG analysis using RNA-sequencing fragments per kilobase of the exon model per million mapped fragment data from the HCC TCGA database. We found a total of 7009 DEGs, of which 5156 genes were up-regulated and 1853 genes were downregulated (Figure 2A). Furthermore, to determine the genes that affect the risk and prognosis of HCC, we performed survival analysis by Cox regression for all 7009 genes. We found that 2667 genes might affect the risk of HCC and its prognosis.

Figure 1

Figure 2

Function analysis of DEGs by Gene Ontology and Kyoto Encyclopedia of Genes and Genomes

To better understand the function of the 2667 DEGs, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis. The GO enrichment analyses. The GO enrichment analysis showed 387 items were related to the DEGs, such as histone−serine phosphorylation and mitotic DNA replication checkpoint. KEGG analysis showed that the DEGs were mainly involved in the nuclear lamina prespliceosome mitotic DNA replication checkpoint. Other pathways included cellular senescence, DNA replication, and oxidation (Figure 2B). Supplementary Tables 2, 3 show more details of GO and KEGG analyses.

Risk scoring model for the application of the LASSO algorithm for HCC

We used LASSO regression to eliminate 2,667 covariates of DEGs that affect the risk of HCC and its prognosis (16). We screened the LASSO regression model using λ = 0.08869376 (Figure 2C). The LASSO regression visualization results are shown in Figure 2D. Finally, 13 hub genes that might predict the prognosis of HCC were obtained. These 13 genes that were DSTNP2, EZH2, GRPEL2, MYCN, NDRG1, PIGU, PPM1G, SF3B4, SQSTM1, TCOF1, TUBG1, TXNRD1, and UCK2. The details of these hub genes are shown in Tables 1, 2. We excluded the pseudogene DSTNP2 because we were more interested in protein expression genes. In addition, we chose eight upregulated genes in HCC and compared with healthy tissues as hub genes. These genes were TXNRD1, TUBG1, SF3B4, PPM1G, PIGU, NDRG1, GRPEL2 and EZH2.

Pan-cancer analysis of DEGs

The eight hub genes upregulated in HCC mentioned above were differentially expressed as shown by pan-cancer analysis (Figure 3A) These eight genes were differentially expressed in digestive tract tumors, such as stomach adenocarcinoma and esophageal carcinoma. GRPEL2 and EZH2 were highly expressed in most cancer types, such as HCC, esophageal carcinoma, and stomach adenocarcinoma. NDRG1 was also highly expressed in HCC, lung squamous cell carcinoma, and head and neck squamous carcinoma. Uterine corpus endometrial carcinoma was associated with genetic mutations in the analyzed samples, and there were only a few mutations in the remaining types of cancer (Figure 3B). In addition, we performed a copy number alteration analysis and showed that most cancer types, inculding HCC, had minor changes in the copy number (Figure 3C). These findings indicated that these eight hub genes were characteristic genes of HCC, and their expression may not be affected by other genomic changes. We then examined the effect of these eight genes on the prognosis of all cancers. We found that these eight genes might affect the prognosis of LGG, HCC, and LUAD. Therefore, a high expression of these genes is a risk factor for these cancers (Figure 3D).

Figure 3

Functional analysis of genes

We conducted a functional analysis of the eight genes described above. First, we used the Gene Set Variation Analysis (GSVA) to calculate the effect of 24 types of immune cells on HCC. This was performed by comparing the correlations between the eight hub genes and immune cells. We found that TUBG1 and SF3B4 could affected most immune cells, and eight of the hub genes were positively related to CD8 T cells and T helper cells. These findings indicated that these eight genes might affect the prognosis of HCC by affecting T helper cells (Figure 4A). Eight of the immune checkpoint genes were related to TNFRSF9, TNFSF9, HAVCR2, and ENTPD1 (Supplementary Table 8 and Supplementary Figure 5). We also analyzed the correlations between the eight hub genes and ferroptosis-related genes. The results showed that the eight genes are positively correlated with most of the ferroptosis-related genes (Supplementary Figure 1, Supplementary Table 5). There was a strong positive correlation between FA complementation group D2(FANCD2), and three characteristic genes, namely EZH2, TUBG1, and PPM1G (Figure 4B). These results indicated that these eight hub genes may affect ferroptosis in HCC. We used a gene set enrichment analysis (GSEA) algorithm to analyze the relationship between the eight hub genes and mitochondria-related pathways. We found that oxidative phosphorylation (OXPHOS), apoptosis, mitophagy, and Fe.S-containing proteins were positively related to these genes according to the enrichment score (Figure 4C). The enrichment scores were positively related to OXPHOS, mitophagy, and Fe.S-containing proteins. On the basis of these results, we hypothesize that ferroptosis affects the prognosis of HCC via mitochondria-related mechanisms.

Figure 4

Validation of gene expression differences and the clinical significance of the prognosis

We measured the expression of the 13 identified hub genes (EZH2, GRPEL2, MYCN, NDRG1, PIGU, PPM1G, DSTNP2, SF3B4, SQSTM1, TCOF1, TUBG1, TXNRD1, and UCK2) by immunohistochemistry (IHC) in tissue samples of 104 patients. We found that eight genes, EZH2, GRPEL2, PIGU, PPM1G, SF3B4, TUBG1, TXNRD1, and NDRG1, were expressed differently between HCC and adjacent non-HCC (Supplementary Table 5, Supplementary Figures 2, 3 and Figure 5A). These results were consistent with bioinformatics analysis. Moreover, the expression levels of these genes were significantly different between HCC and adjacent non-HCC lesions (Supplementary Table 5). In conclusion, these eight genes could be used as potential risk factor indicators to determine clinical indications for prognosis. To further ensure the accuracy of core genes, we used TCGA+GTEx and GSE101685 data to further analyze the differential expression of these genes. We found that the eight hub genes identified above were all differentially expressed (Figure 6 and Supplementary Table 7).

Figure 5

Figure 6

HCC prognosis-related model construction

In addition to analyzing expression of the eight hub genes, we collected basic information on all 104 patients, such as sex, age, histological differentiation, vascular cancer embolus, cirrhosis, tumor size (cm), alpha-fetoprotein concentrations (ng/mL), hepatitis B surface antigen concentrations, alanine aminotransferase, underlying diseases, and total bilirubin concentrations. We excluded hepatitis C viruse infection patients in our study. These clinical factors are generally considered indicators of the prognosis of HCC. In our dataset, the one-factor survival analysis showed that a larger tumor size and positive vascular cancer embolus indicated a poor prognosis. Therefore, we incorporated these factors into the Cox regression model in the multivariate analysis. The high expression of EZH2, GRPEL2, TCOF1 and NDRG1 genes suggested a poor prognosis, as shown in the multivariate regression analysis (Table 3, Figure 5B). Furthermore, a nomogram was used for visualization to better use as a clinical prognostic model. The nomogram showed that a larger tumor size and higher EZH2, GRPEL2, and NDRG1expression predicted the higher likelihood of a worse prognosis (Figure 5C). The nomogram was able to predict the probability of 1- and 3-year survival (Figure 5D). These findings indicated a close relation to survival, which suggested that our nomogram may be a reliable prognostic model for HCC. Furthermore, we computed the C-index to evaluate our model. The C-index of this model was 0.811(0.797-0.825). A decision curve showed that this model had reliable clinical practicability (Supplementary Figure 4).

HCC prognosis-related model functional analysis

To further understand the function of the HCC prognostic model, we analyzed the relationship between model scores and clinical parameters. We found that the prognostic model score was mainly related to alpha-fetoprotein, T stage. Among them, the model score is higher when the disease is more severe (Figure 7A). In addition, we analyzed the relationship between 3 genes and HBV infection. EZH2 was related to HBV infection, and the expression of EZH2 expression was high in patients with HBV infection (Supplementary Table 9). In addition, GSVA scores were calculated for tumor-related pathways, and the correlations between model scores and tumor-related pathways were determined. We found that the HCC model was mainly related to the inflammatory response and other pathways (Figure 7B).

Figure 7

Prognostic models and drug resistance

Our analysis using GSE106988 showed that NDRG1 was mainly associated with functions such as RNA shearing(Figure 8A), but the other two genes need to be verified in subsequent experiments. After analyzing the targets of anti-tumor drug sensitivity action in patients with HCC, we found that the prognostic model of ferroptosis was associated with multiple chemotherapy efficacy, with the strongest correlation with axitinib (Figure 8B). We also analyzed the function of the three core genes for ferroptosis and the regulatory relationships. Finally, to understand the post-transcriptional regulatory relationships of these three genes, analysis using POSTAR3 analysis showed that these three genes were affected by multiple post-transcriptional regulatory genes, among which IGF2BP2, FIPILI, TARDBP, ELF2, ATXN2, DDx3X, AGO2, IGF2BP1, and HNRNPC regulated these three genes simultaneously (Figure 8C).

Figure 8

Discussion

HCC is a common and highly malignant tumor. Patients with HCC have a poor prognosis (17). At present, surgery is the most effective treatment method for HCC (18). However, at diagnosis, most patients with HCC show intra- and extra-hepatic metastases, vascular tumor thrombi and inadequate future liver volume, which indicate that the opportunity for surgery has already passed. With the development of advanced therapies, various emerging treatments have achieved positive results (19, 20). However, because of the complexity of the genetic background of HCC, the benefits of these treatments remain unclear (20, 21). Therefore, in-depth studies of the pathogenesis of HCC, including machine learning, will have clinical significance for predicting the interaction of hub genes and investgating prognostic-related molecular markers for individuals. In this study, we constructed a prognostic model combined with clinical information, functional analysis, and pan-cancer analyses to verify our bioinformatic results (22).

In this study, we used the TCGA data sets to predict cancer the risk of HCC and its prognosis. We found that the DEGs were highly related to DNA replication, oxidation, and inflammatory response. According to previous studies, these mechanisms are consistent with the activity of tumor cell growth and mitochondrial metabolism (23). Therefore, we speculate that abnormal mitochondrial function is related to the occurrence and development of HCC. The LASSO regression analysis was performed to generate a prognostic signature. We found 13 hub genes that were potentially related to the prognosis of HCC. Eight upregulated genes (EZH2, GRPEL2, PIGU, PPM1G, SF3B4, TUBG1, TXNRD1, and NDRG1) were not only found to be tumorigenesis risk factors, but also factors of a poor prognosis of HCC. In addition, we analyzed 104 clinical samples to validate our predictions. These eight genes were differentially expressed between HCC and adjacent non-HCC lesions, which suggested that these genes could be used as biomarkers for risk factors. The multivariate regression analysis showed that EZH2, NDRG1 and GRPEL2, as well as the tumor size were independent factors for the prognosis of HCC. According to recent surveys, TOCF1 (11) and UCK2 (24) may also indicate the prognosis in the patients with HCC. Unfortunately, our study did not show any statistical significance for the differentiated expression of TOCF1 or UCK2. Therefore, a larger sample size and multicenter study would be useful to further validate the relationship between the expression of these two genes and the prognosis of HCC.

We also constructed a prognostic nomogram model, which included EZH2, GRPEL2, NDRG1, and the tumor size. EZH2 is a transcriptional suppressor gene (25). EZH2 combines with EED, SUZ12, and RBBP4 to generate a PRC2 subunit, which is overexpressed in multiple types of cancer (26). EZH2 has a bidirectional tumorigenic effect involving mTOR signaling, as shown in EZH2-deficient leukemia cells (27). In addition, detecting EZ gene expression by tissue puncture biopsy can be used as a molecular marker for the prognosis of patients with HCC (28). There have been several recent clinical trials of EZH2 inhibitors in various cancers being studied (29–32). Our study showed that EZH2 was involved in autophagy, which is consistent with previous studies. EZH2 may also have specific effects on mitophagy, ferroptosis, and OXPHOS. Determining the specific mechanism of EZ2H should be considered in future research. NDRG1 plays a crucial role in angiogenesis, proliferation, invasion, and prognosis of HCC (33, 34). NDRG1 may also be a molecular marker for metastasis and prognosis in HCC. NDRG1 may predict the recurrence of HCC in liver transplant patients after surgery (35). This gene plays an anti-tumor role by involving the Wnt signaling pathway and other carcinogenesis signaling pathways. We speculate that NDRG1 may be involved in metabolic processes related to ferroptosis. However, more experimental studies still need to verify the mechanism of NDRG1 in ferroptosis. Overexpression of GRPEL2 affects the stability of mitochondrial proteins, and GRPEL2 may protect protein stability under mitochondrial oxidative stress (36). GRPEL2 is responsible for maintaining mitochondrial homeostasis through the nuclear factor-kappa B pathway, which affects the cell cycle in HCC (37). In particular, the composition of Fe.S-containing proteins is related to DEGs, may have an anti-tumor effect via mitochondrial-related pathways, such as complex I, OXPHOS, and mitophagy. Our finding suggested that upregulated NDRG1 could be related to poor prognosis through mitochondrial-related proteins in HCC. On the basis of the results of available datasets, we found that NDRG1 silencing was associated with various tumor biological functions such as RNA splicing, histone modification, and peptidyl-lysine modification. We concluded that the ferroptosis-related model was associated with the prognosis of HCC. These 3 genes may be essential targets and regulatory proteins in the biological progression of HCC. Using drug sensitivity analysis, we found a significant negative association of multiple anti-cancer drugs with the model. Multiple post-transcriptional regulatory mechanisms may be involved. Among them, Axitinib is an anti-tumor agent used clinically in adult patients with progressive renal cell carcinoma (RCC) who have failed prior treatment with a tyrosine kinase inhibitor or cytokine. This study provides some theoretical reference to determine whether Axitinib can be administered to patients with HCC based on ferroptosis-related prediction models. However clinical trials are still required to validate this theory. We also included a clinical factor, tumor size, in the prognosis model. Patients with more extensive tumor volumes have a worse prognosis, because the tumor has adapted to the body’s growth or because the tumor is discovered late. Therefore, the current research results suggest that EZH2, GRPEL2, NDRG1, and tumor size are reliable prognostic indicators.

The pan-cancer analysis results suggest that the eight hub genes (TXNRD1, TUBG1, SF3B4, PPM1G, PIGU, NDRG1, GRPEL2, and EZH2) were differentially expressed in gastrointestinal tumors. These eight genes may have unique roles in the digestive tract and are relatively stable in genomics. Cirrhosis-related HCC characteristics of HBV infection have a particular tumor immune microenvironment. This microenvironment may be in a state of immunosuppression over a long period, may play a unique role in promoting HCC to escape host immune surveillance (38). Therefore, we speculate that prognosis-related genes could affect the prognosis through the tumor immune microenvironment in HCC. However, further studies are required to determine the mechanisms of interaction between hub genes and immune cells. We also wanted to determine the role of the tumor microenvironment and energy metabolism in HCC. One of our study aim was to detect an unknown mechanism that links programmed cell death, and mitochondrial metabolism in HCC. Therefore, we performed an immune infiltration analysis and a GSEA analysis. In the immune infiltration analysis, we found that expression of the eight hub genes (TXNRD1, TUBG1, SF3B4, PPM1G, PIGU, NDRG1, GRPEL2, EZH2) and the immune infiltration of T helper cells, Tcm cells, and Th2 cells were related to HCC. Studies have shown that the infiltration of multiple immune cells plays a role in HBV-related HCC (39). T cells may mediate cell death by various mechanisms. Previous studies have demonstrated that Th17 cells may activate and induce cell death by Fas-mediated apoptosis, which is consistent with our results (40, 41).

Ferroptosis is an iron-dependent regulatory form of cell death caused by excessive lipid peroxidation, which is associated with the occurrence and response to treatment of various types of tumors. Ferroptotic injury triggers inflammation-related immunosuppression in the tumor microenvironment, which favors tumor growth. Ferroptosis can also lead to the activation of the Renin-angiotensin system, TP53 and other signaling pathways, which then leads to poor survival of patients with HCC (42). We constructed a model in which EZH2, TUBG1, and PPM1G were strongly correlated with the FANCD2. FANCD2 is an intersection protein associated with multiple ferroptosis in clear cell renal cell carcinoma and bladder cancer (43, 44). According to our results, EZH2, TUBG1, and PPM1G may play a role in ferroptosis through the key gene FANCD2. Fe.S-containing proteins were also significantly correlated with the eight hub genes. Current research has shown that Fe.S proteins participate in vital pathways such as oxidative phosphorylation and iron metabolism. Fe.S proteins are involved in the critical steps of the pathway that involve the maturation of mitochondrial proteins. The vital predisposing factors for ferroptosis are iron leakage, lipid peroxidation, and increased ROS (45, 46). The suppression of ferroptosis may contribute to tumor progression and survival (47). Some scholars believe that mitochondria play a crucial role in ferroptosis caused by erastin or cystine deprivation. However, a lack of mitochondria is affected by erastin concentrations, leading to cell death, which suggests that mitochondria may participate in ferroptosis through other currently unknown mechanisms (48). Because of that, ferroptosis mediates various tumor effects. In different damage-associated molecular patterns, ferroptosis can protect cells from drug damage (42). Such mechanisms may aggravate the damage caused by immune cells (49). Therefore, we speculate that the development of HCC may be affected by ferroptosis via mitochondrial metabolism and ROS. Additional mechanistic studies are still required to confirmed this possibility.

On the basis of the above-mentioned verification results of our model’s function, ferroptosis, mitochondrial metabolism, ROS, and immune infiltration are closely related to the inflammatory background of HCC. The background of hepatitis B infection and chronic inflammation exhibits a unique immune microenvironment (50, 51). Inflammation is tightly associated with cell death (52). Functional studies have shown that clinical parameters and inflammation are positively correlated with the prognostic models. Programmed cell death is affected by chronic inflammation, which results in tissue DNA damage (53). This study showed that our model was positively correlated with inflammation and DNA damage repair. Therefore, this model can reflect the inflammation associated with cell death in HCC. mTOR complex I can regulate endoplasmic reticulum (ER) stress (54). Transcription factor families, such as Myc and E2F, target various tumorigenesis and mediate inflammatory progression (55). Therefore, our model can reflect the prognosis of HCC in the context of inflammation.

There are some limitations to our study. First, this study recruited 104 patients with HCC in East China. This sample size was limited to a specific race and a single center. To apply our findings in the clinic, multicenter studies with a larger sample size are required. In particular, whether our prognostic model can be used as a biomarker in a large sample size needs to be investigated. Second, we developed a prognostic model for cell death associated with ferroptosis. However, several mechanisms have not yet been fully determined, and more evidence is required. In particular, ferroptosis-related therapeutic regimens should be designed for these targets and validated in in vivo and in vitro studies. Finally, further experiments are still required to verify the mechanism of action of FANCD2 through ferroptosis affecting the inflammatory and immune microenvironment of HCC. The most appropriate treatment required for patients with a high risk of recurrence and the most suitable types of surgery or adjuvant treatment required to prolong patients’ overall survival are still unknown.

In conclusion, we used bioinformatics analysis methods to predict DEGs for the occurrence and development of HCC and found eight hub genes (EZH2, GRPEL2, PIGU, PPM1G, SF3B4, TUBG1, TXNRD1, and NDRG1). We further verified the clinical tissue sample results to develop a reliable nomogram, which could be used as a prognostic model for HCC. These hub genes are involved in mitophagy, OXPHOS, and Fe.S-containing proteins. EZH2, TUBG1, and PPM1G are related to the ferroptosis gene, FANCD2, and may be involved in tumor biological behaviors mediated by ferroptosis. The background of chronic hepatitis may play an essential role in the development of HCC. According to our data set, three genes, EZH2, GRPEL2, and NDRG1, and one clinical factor, tumor size, could be used as reliable indicators for evaluating the prognosis in patients with HCC.

References

• 1

  EstesCRazaviHLoombaRYounossiZSanyalAJ. Modeling the epidemic of nonalcoholic fatty liver disease demonstrates an exponential increase in burden of disease. Hepatology (2018) 67:123–33. doi: 10.1002/hep.29466

  • CrossRef
  • Google Scholar

• 2

  MichielsenPPFrancqueSMvan DongenJL. Viral hepatitis and hepatocellular carcinoma. World J Surg Oncol (2005) 3:27. doi: 10.1186/1477-7819-3-27

  • CrossRef
  • Google Scholar

• 3

  BrayFFerlayJSoerjomataramISiegelRLTorreLAJemalA. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2018) 68:394–424. doi: 10.3322/caac.21492

  • CrossRef
  • Google Scholar

• 4

  WuJLuWYCuiLL. Clinical significance of STAT3 and MAPK phosphorylation, and the protein expression of cyclin D1 in skin squamous cell carcinoma tissues. Mol Med Rep (2015) 12:8129–34. doi: 10.3892/mmr.2015.4460

  • CrossRef
  • Google Scholar

• 5

  WangLWuJSongSChenHHuYXuBet al. Plasma exosome-derived sentrin SUMO-specific protease 1: A prognostic biomarker in patients with osteosarcoma. Front Oncol (2021) 11:625109. doi: 10.3389/fonc.2021.625109

  • CrossRef
  • Google Scholar

• 6

  RizzolloFMoreSVangheluwePAgostinisP. The lysosome as a master regulator of iron metabolism. Trends Biochem Sci (2021) 46(12):960–75. doi: 10.1016/j.tibs.2021.07.003

  • CrossRef
  • Google Scholar

• 7

  YeZLiuWZhuoQHuQLiuMSunQet al. Ferroptosis: Final destination for cancer? Cell Prolif (2020) 53:e12761. doi: 10.1111/cpr.12761

  • CrossRef
  • Google Scholar

• 8

  LiuYGuW. The complexity of p53-mediated metabolic regulation in tumor suppression. Semin Cancer Biol (2021) 85:4–32. doi: 10.1016/j.semcancer.2021.03.010

  • CrossRef
  • Google Scholar

• 9

  SuYZhaoBZhouLZhangZShenYLvHet al. Ferroptosis, a novel pharmacological mechanism of anti-cancer drugs. Cancer Lett (2020) 483:127–36. doi: 10.1016/j.canlet.2020.02.015

  • CrossRef
  • Google Scholar

• 10

  WanSLeiYLiMWuB. A prognostic model for hepatocellular carcinoma patients based on signature ferroptosis-related genes. Hepatol Int (2021) 16(1):112–24. doi: 10.21203/rs.3.rs-614871/v1

  • CrossRef
  • Google Scholar

• 11

  WangZFuYXiaAChenCQuJXuGet al. Prognostic and predictive role of a metabolic rate-limiting enzyme signature in hepatocellular carcinoma. Cell Prolif (2021) 54(10):e13117. doi: 10.1111/cpr.13117

  • CrossRef
  • Google Scholar

• 12

  LiPZhangYXuYCaoHLiLXiaoH. Characteristics of CD8+ and CD4+ tissue-resident memory lymphocytes in the gastrointestinal tract. Advanced Gut Microbiome Res (2022) 2022:1–12. doi: 10.1155/2022/9157455

  • CrossRef
  • Google Scholar

• 13

  Gene OntologyC. The gene ontology (GO) project in 2006. Nucleic Acids Res 34 (2006) 34(Database issue):D322–6. doi: 10.1093/nar/gkj021

  • CrossRef
  • Google Scholar

• 14

  KanehisaMGotoSKawashimaSOkunoYHattoriM. The KEGG resource for deciphering the genome. Nucleic Acids Res (2004) 32:D277–80. doi: 10.1093/nar/gkh063

  • CrossRef
  • Google Scholar

• 15

  TibshiraniR. The lasso method for variable selection in the cox model. Stat Med (1997) 16:385–95. doi: 10.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3

  • CrossRef
  • Google Scholar

• 16

  SimonNFriedmanJHastieTTibshiraniR. Regularization paths for cox's proportional hazards model via coordinate descent. J Stat Softw (2011) 39:1–13. doi: 10.18637/jss.v039.i05

  • CrossRef
  • Google Scholar

• 17

  WenNCaiYLiFYeHTangWSongPet al. The clinical management of hepatocellular carcinoma worldwide: A concise review and comparison of current guidelines: 2022 update. Biosci Trends (2022) 16:20–30. doi: 10.5582/bst.2022.01061

  • CrossRef
  • Google Scholar

• 18

  MidorikawaYTakayamaTHigakiTAramakiOTeramotoKYoshidaNet al. Comparison of the surgical outcomes in patients with synchronous versus metachronous multiple hepatocellular carcinoma. Biosci Trends (2021) 14:415–21. doi: 10.5582/bst.2020.03313

  • CrossRef
  • Google Scholar

• 19

  HaradaMAramakiOMidorikawaYHigakiTNakayamaHMoriguchiMet al. Impact of patient age on outcome after resection for hepatocellular carcinoma. Biosci Trends (2021) 15:33–40. doi: 10.5582/bst.2020.03437

  • CrossRef
  • Google Scholar

• 20

  LiuDSongT. Changes in and challenges regarding the surgical treatment of hepatocellular carcinoma in China. Biosci Trends (2021) 15:142–7. doi: 10.5582/bst.2021.01083

  • CrossRef
  • Google Scholar

• 21

  TangHCaoYJianYLiXLiJZhangWet al. Conversion therapy with an immune checkpoint inhibitor and an antiangiogenic drug for advanced hepatocellular carcinoma: A review. Biosci Trends (2022) 16:130–41. doi: 10.5582/bst.2022.01019

  • CrossRef
  • Google Scholar

• 22

  WuJChenZPShangAQWangWWChenZNTaoYJet al. Systemic bioinformatics analysis of recurrent aphthous stomatitis gene expression profiles. Oncotarget (2017) 8:111064–72. doi: 10.18632/oncotarget.22347

  • CrossRef
  • Google Scholar

• 23

  XiaoZDaiZLocasaleJW. Metabolic landscape of the tumor microenvironment at single cell resolution. Nat Commun (2019) 10:3763. doi: 10.1038/s41467-019-11738-0

  • CrossRef
  • Google Scholar

• 24

  LiYQiDZhuBYeX. Analysis of m6A RNA methylation-related genes in liver hepatocellular carcinoma and their correlation with survival. Int J Mol Sci (2021) 22(3):1474. doi: 10.3390/ijms22031474

  • CrossRef
  • Google Scholar

• 25

  GuoJDaiXLaurentBZhengNGanWZhangJet al. AKT methylation by SETDB1 promotes AKT kinase activity and oncogenic functions. Nat Cell Biol (2019) 21:226–37. doi: 10.1038/s41556-018-0261-6

  • CrossRef
  • Google Scholar

• 26

  MargueronRReinbergD. The polycomb complex PRC2 and its mark in life. Nature (2011) 469:343–9. doi: 10.1038/nature09784

  • CrossRef
  • Google Scholar

• 27

  GuZLiuYCaiFPatrickMZmajkovicJCaoHet al. Loss of EZH2 reprograms BCAA metabolism to drive leukemic transformation. Cancer Discov (2019) 9:1228–47. doi: 10.1158/2159-8290.CD-19-0152

  • CrossRef
  • Google Scholar

• 28

  CaiMYTongZTZhengFLiaoYJWangYRaoHLet al. EZH2 protein: a promising immunomarker for the detection of hepatocellular carcinomas in liver needle biopsies. Gut (2011) 60:967–76. doi: 10.1136/gut.2010.231993

  • CrossRef
  • Google Scholar

• 29

  YinHWangYWuYZhangXZhangXLiuJet al. EZH2-mediated epigenetic silencing of miR-29/miR-30 targets LOXL4 and contributes to tumorigenesis, metastasis, and immune microenvironment remodeling in breast cancer. Theranostics (2020) 10:8494–512. doi: 10.7150/thno.44849

  • CrossRef
  • Google Scholar

• 30

  RizqOMimuraNOshimaMSarayaAKoideSKatoYet al. Dual inhibition of EZH2 and EZH1 sensitizes PRC2-dependent tumors to proteasome inhibition. Clin Cancer Res (2017) 23:4817–30. doi: 10.1158/1078-0432.CCR-16-2735

  • CrossRef
  • Google Scholar

• 31

  LangJDHendricksWPDOrlandoKAYinHKieferJRamosPet al. Ponatinib shows potent antitumor activity in small cell carcinoma of the ovary hypercalcemic type (SCCOHT) through multikinase inhibition. Clin Cancer Res (2018) 24:1932–43. doi: 10.1158/1078-0432.CCR-17-1928

  • CrossRef
  • Google Scholar

• 32

  ZhangPXiaoZWangSZhangMWeiYHangQet al. ZRANB1 is an EZH2 deubiquitinase and a potential therapeutic target in breast cancer. Cell Rep (2018) 23:823–37. doi: 10.1016/j.celrep.2018.03.078

  • CrossRef
  • Google Scholar

• 33

  NishioSUshijimaKTsudaNTakemotoSKawanoKYamaguchiTet al. Cap43/NDRG1/Drg-1 is a molecular target for angiogenesis and a prognostic indicator in cervical adenocarcinoma. Cancer Lett (2008) 264:36–43. doi: 10.1016/j.canlet.2008.01.020

  • CrossRef
  • Google Scholar

• 34

  YanXChuaMSSunHSoS. N-myc down-regulated gene 1 mediates proliferation, invasion, and apoptosis of hepatocellular carcinoma cells. Cancer Lett (2008) 262:133–42. doi: 10.1016/j.canlet.2007.12.010

  • CrossRef
  • Google Scholar

• 35

  ChengJXieHYXuXWuJWeiXSuRet al. NDRG1 as a biomarker for metastasis, recurrence and of poor prognosis in hepatocellular carcinoma. Cancer Lett (2011) 310:35–45. doi: 10.1016/j.canlet.2011.06.001

  • CrossRef
  • Google Scholar

• 36

  KonovalovaSLiuXManjunathPBaralSNeupaneNHilanderTet al. Redox regulation of GRPEL2 nucleotide exchange factor for mitochondrial HSP70 chaperone. Redox Biol (2018) 19:37–45. doi: 10.1016/j.redox.2018.07.024

  • CrossRef
  • Google Scholar

• 37

  LaiMCZhuQQXuJZhangWJ. Experimental and clinical evidence suggests that GRPEL2 plays an oncogenic role in HCC development. Am J Cancer Res (2021) 11(9):4175–98.

  • Google Scholar

• 38

  BonnelARBunchorntavakulCReddyKR. Immune dysfunction and infections in patients with cirrhosis. Clin Gastroenterol Hepatol (2011) 9:727–38. doi: 10.1016/j.cgh.2011.02.031

  • CrossRef
  • Google Scholar

• 39

  OuraKMorishitaATaniJMasakiT. Tumor immune microenvironment and immunosuppressive therapy in hepatocellular carcinoma: A review. Int J Mol Sci 22 (2021) 22(11):5801. doi: 10.3390/ijms22115801

  • CrossRef
  • Google Scholar

• 40

  YuYIclozanCYamazakiTYangXAnasettiCDongCet al. Abundant c-fas-associated death domain-like interleukin-1-converting enzyme inhibitory protein expression determines resistance of T helper 17 cells to activation-induced cell death. Blood (2009) 114:1026–8. doi: 10.1182/blood-2009-03-210153

  • CrossRef
  • Google Scholar

• 41

  BortnerCDScoltockABCainDWCidlowskiJA. T-Cell development of resistance to apoptosis is driven by a metabolic shift in carbon source and altered activation of death pathways. Cell Death Differ (2016) 23:889–902. doi: 10.1038/cdd.2015.156

  • CrossRef
  • Google Scholar

• 42

  ChenXKangRKroemerGTangD. Broadening horizons: the role of ferroptosis in cancer. Nat Rev Clin Oncol (2021) 18:280–96. doi: 10.1038/s41571-020-00462-0

  • CrossRef
  • Google Scholar

• 43

  WuGWangQXuYLiQChengL. A new survival model based on ferroptosis-related genes for prognostic prediction in clear cell renal cell carcinoma. Aging (Albany NY) (2020) 12:14933–48. doi: 10.18632/aging.103553

  • CrossRef
  • Google Scholar

• 44

  LiuJMaHMengLLiuXLvZZhangYet al. Construction and external validation of a ferroptosis-related gene signature of predictive value for the overall survival in bladder cancer. Front Mol Biosci (2021) 8:675651. doi: 10.3389/fmolb.2021.675651

  • CrossRef
  • Google Scholar

• 45

  WangYWeiZPanKLiJChenQ. The function and mechanism of ferroptosis in cancer. Apoptosis (2020) 25:786–98. doi: 10.1007/s10495-020-01638-w

  • CrossRef
  • Google Scholar

• 46

  DixonSJLembergKMLamprechtMRSkoutaRZaitsevEMGleasonCEet al. Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell (2012) 149:1060–72. doi: 10.1016/j.cell.2012.03.042

  • CrossRef
  • Google Scholar

• 47

  XieYHouWSongXYuYHuangJSunXet al. Ferroptosis: process and function. Cell Death Differ (2016) 23:369–79. doi: 10.1038/cdd.2015.158

  • CrossRef
  • Google Scholar

• 48

  ZhengJConradM. The metabolic underpinnings of ferroptosis. Cell Metab (2020) 32:920–37. doi: 10.1016/j.cmet.2020.10.011

  • CrossRef
  • Google Scholar

• 49

  ZdralevicMVuceticMDaherBMarchiqIParksSKPouyssegurJ. Disrupting the 'Warburg effect' re-routes cancer cells to OXPHOS offering a vulnerability point via 'ferroptosis'-induced cell death. Adv Biol Regul (2018) 68:55–63. doi: 10.1016/j.jbior.2017.12.002

  • CrossRef
  • Google Scholar

• 50

  XiaYCarpentierAChengXBlockPDZhaoYZhangZet al. Human stem cell-derived hepatocytes as a model for hepatitis b virus infection, spreading and virus-host interactions. J Hepatol (2017) 66:494–503. doi: 10.1016/j.jhep.2016.10.009

  • CrossRef
  • Google Scholar

• 51

  KhanamAChuaJVKottililS. Immunopathology of chronic hepatitis b infection: Role of innate and adaptive immune response in disease progression. Int J Mol Sci (2021) 22(11):5497. doi: 10.3390/ijms22115497

  • CrossRef
  • Google Scholar

• 52

  CiccaroneFCastelliSCirioloMR. Oxidative stress-driven autophagy acROSs onset and therapeutic outcome in hepatocellular carcinoma. Oxid Med Cell Longev (2019) 2019:6050123. doi: 10.1155/2019/6050123

  • CrossRef
  • Google Scholar

• 53

  Malcov-BrogHAlpertAGolanTParikhSNordlingerANettiFet al. UV-Protection timer controls linkage between stress and pigmentation skin protection systems. Mol Cell (2018) 72:444–456 e7. doi: 10.1016/j.molcel.2018.09.022

  • CrossRef
  • Google Scholar

• 54

  Appenzeller-HerzogCHallMN. Bidirectional crosstalk between endoplasmic reticulum stress and mTOR signaling. Trends Cell Biol (2012) 22:274–82. doi: 10.1016/j.tcb.2012.02.006

  • CrossRef
  • Google Scholar

• 55

  WenZJinKShenYYangZLiYWuBet al. N-myristoyltransferase deficiency impairs activation of kinase AMPK and promotes synovial tissue inflammation. Nat Immunol (2019) 20:313–25. doi: 10.1038/s41590-018-0296-7

  • CrossRef
  • Google Scholar