The gene expression microarray and corresponding clinical information were acquired from The Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC, https://cancergenome.nih.gov/). To validate the prognostic potential of the genetic risk score, the independent HCC datasets were obtained through the International Cancer Genome Consortium (ICGC, https://icgc.org) and the Gene Expression Omnibus (GEO) database (GSE14520, Affymetrix HT Human Genome U133A Array, GPL3921). The protein expression profiles of HCC were acquired from The National Cancer Institute’s Clinical Proteomic Tumor Analysis Consortium (CPTAC, https://cptac-data-portal.georgetown.edu/). The ferroptosis-related genes (FRGs) were acquired from the FerrDb website (http://www.zhounan.org/ferrdb). All statistical analyses were performed in R version 4.0.5.

Differentially expressed genes (DEGs) were determined using R package “limma” (Ritchie et al., 2015) based on the comparison of gene expression levels between HCC and liver tissues in TCGA cohort. The selection criteria for identification of DEGs were as follows: |logFC|>1 and adjusted p < 0.01. This study fully complies with TCGA publication requirements (http://cancergenome.nih.gov/publications/publicationguidelines). Kaplan–Meier survival analysis and the Cox proportional hazards model were used to analyze the association between the different clusters and prognosis with the R package “Survminer.” The receiver operating characteristic (ROC) curve was used to assess the prognosis classification performance of the signature and nomogram, and the area under the curve (AUC) was calculated using R package “timeROC.”

The R package “clusterProfiler” was used for Gene Ontology (GO) and KEGG pathway analysis (Yu et al., 2012). p < 0.01 and FDR <0.05 were considered statistically significant. Gene set enrichment analysis (GSEA) of hallmark gene sets was acquired from MSigDB database v7.1 (Hänzelmann et al., 2013; Subramanian et al., 2005).

The genomic mutation data were curated from TCGA database. The copy number variation distribution information of FRGs on different chromosomes was plotted using the R package “Rcircos.”

To study the function of DEGs, we separated patients into three subclusters by the R package of “ConsensusClusterPlus” (Wilkerson and Hayes, 2010). The R package “ggplot2” of principal component analysis (PCA) was carried out to explore the genes in subclusters of HCC.

The enrichment scores of immune infiltration levels were calculated by single sample gene set enrichment analysis (ssGSEA). ssGSEA is a popular enrichment algorithm, which was extensively utilized in medical studies (Liu et al., 2021a; Liu et al., 2021b; Liu et al., 2022a; Liu et al., 2022b; Liu et al., 2022c). Data were visualized using multidimensional scaling (MDS) and a Gaussian fitting model. The immune infiltration assessment was performed using the Tumor IMmune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/) and microenvironment cell population count (MCP-counter) method.

In the present study, a nomogram was constructed based on the risk model and pathological stage to assess the probability of individual patient’s OS with HCC. Meanwhile, the calibration curves were performed to verify the accuracy of the nomogram. Furthermore, the ROC and DCA were applied to evaluate the prediction efficiency and net benefit of the nomogram for HCC patients at 1, 3, and 5 years.

A workflow of this study is presented in Figure 1. Based on the threshold of FDR <0.05 and |log2 FC| > 1, a total of 51 upregulated and 7 downregulated significant ferroptosis-related differentially expressed genes (DEGs) were identified (Figure 2A), and the expression landscape of DEGs was shown in the heatmap (Figure 2B). Next, the ferroptosis-related DEGs were selected to perform the GO (Figure 2C) and KEGG (Figure 2D) analyses. Subsequently, several biological processes, including response to oxidative stress, reactive oxygen species metabolic process, and aging, were detected in GO terms. As to KEGG pathway analysis, ferroptosis, serotonergic synapse, VEGF signaling pathway, arachidonic acid metabolism, etc. were mostly associated with these genes.

Univariate Cox regression analysis was performed to evaluate the prognostic value of ferroptosis-related DEGs, and 25 DEG aberrant expressions were correlated with the overall survival of HCC patients in TCGA cohort (Figure 3A). Mutation analysis of 25 prognostic genes revealed that 10% of genes were mutated by cBioportal (Figure 3B). We first figured out the somatic mutation rates of 25 DEGs based on the simple nucleotide variation data from TCGA database. The results revealed that only 39 of the 364 (10.71%) samples presented genetic alterations. We observed CNV changed in YY1AP1, STMN1, CDKN2A, and SLC2A4 (Figure 3C). The location of 25 ferroptosis-related genes on the chromosomes is shown in Figure 3D. The PPI network of ferroptosis-related DEGs is constructed in Figure 3E.

The “ConcensusCluster” package in R language was used to perform consensus clustering analyses of HCC samples from TCGA dataset based on the 25 ferroptosis-related genes shown in Figures 4A,B. In order to obtain a reliable and stability cluster, k = 3 was chosen from k = 2 to 9 in TCGA dataset. However, the PCA analysis revealed that k = 3 was also preferable (Figure 4C). Then, the HCC patients were split into three distinct clusters: A, B, and C clusters. Figure 4D shows an obvious distinction in overall survival among these three clusters. K–M analysis suggested that significant overall survival benefit of the A cluster over B or C cluster, and B cluster over C cluster (Figures 4D–G). After combining clusters A and B, a significantly longer overall survival was shown for A + B than for C by log-rank tests (p < 0.001; Figure 4H). In addition, after combining clusters B and C, a significantly longer overall survival was shown for A than for B + C by log-rank tests (p < 0.001; Figure 4I). A clear disease-free survival (DFS) difference among the three clusters was observed (p < 0.001; Figure 5A). The DFS was longer for cluster A than for cluster B patients, and cluster B patients than for cluster C patients. The DFS in cluster A was significantly longer than that in clusters B and C (p < 0.001, respectively; Figure 5B,C). The DFS was longer for cluster B than for cluster C patients (p < 0.001; Figure 5D). After combining clusters A and B, a significantly longer overall survival was shown for A + B than for C by log-rank tests (p < 0.001; Figure 5E). In addition, after combining clusters B and C, a significantly longer overall survival was shown for A than for B + C by log-rank tests (p < 0.001; Figure 5F). Figures 5G,H indicate that the three clusters had different clinicopathologic features, including T stage, tumor stage, grade, gender, age, and survival status.

Moreover, the ICGC cohort was regarded as the external validation set. Figure 6A shows an obvious distinction in overall survival among these three clusters (p = 0.022). Figures 6B,C show patients in cluster A had longer OS time than patients in cluster B (p = 0.009) and cluster C (p = 0.019). But the survival difference between clusters B and C was not significant (p = 0.912) (Figure 6D). After combining clusters B and C, patients in cluster A presented longer OS time than patients in clusters B + C (Figure 6E; p = 0.006). But the survival difference between clusters A + B and C was not significant (Figure 6F; p = 0.275). The analysis further showed there were significant differences in the TNM stage and survival status (Figures 6G,H).

Immune cells are well known to be the key part for antitumor immune response and immunotherapy; therefore, we examined the distribution of immune cell infiltration in the three clusters based on the ssGSEA algorithm (Figure 7A). A higher ssGSEA score indicated more infiltrating immune cells. The results showed that activated CD4+/CD8+ T cells, effector memory CD4+/CD8+ T cells, activated dendritic cells, activated B cells, immature B cells, immature dendritic cells, MDSC, regulatory T cells, T follicular helper cells, and type 17 T helper cells were mainly enriched in the C subtype. However, eosinophils were markedly enriched in the A subtype. Furthermore, we also assessed the immune infiltration level using TIMER and MCPcount algorithms, which indicated a significant difference between the three clusters (Figures 7B,C; p < 0.05). In addition, there was clear difference in the expression of PD-L1 and TP53 between the three clusters (Figures 7D,E; p < 0.05).

To explore the biological molecular changes among these ferroptosis-related clusters, we performed GSVA enrichment. Heatmap showed that P53 pathway, DNA replication, and cell cycle were more enriched in cluster B than in cluster A (Figure 8A); cell cycle, DNA replication, and mismatch repair were more enriched in cluster C than in cluster A (Figure 8B); the and PPAR pathway, fatty acid metabolism, and multiple metabolism pathway were more enriched in cluster C than in cluster B (Figure 8C). The color changes from yellow to red, determining an increase in the value of the enriched score. Blue represents cluster A, yellow represents cluster B, and red represents cluster C. GSVA showed that multiple metabolic and cell cycle pathways were identified. A total of 1,966 differential genes were obtained among the three clusters (Figure 8D). Functional enrichment analysis was performed through the GO and KEGG (Figures 8E,F).

Based on the univariate Cox regression analysis of TCGA cohort, genes with significant prognostic value were further extracted (p < 0.01; Figure 9A). By the random survival forest variable hunting (RSFVH) algorithm, the top 10 genes were ranked based on their predictive importance: TMEM251, MED8, UCK2, STIP1, PIGU, PPM1G, ZDHHC18, MUTYH, RAN, and SNRPB (Figure 9B). We performed Kaplan–Meier analysis to figure out the best risk model with relatively small genes and relative significant p value. Then, we found the five-gene signature ranked top, including MED8, PIGU, PPM1G, RAN, and SNRPB (Figure 9C). Subsequently, we inspected the distribution of risk score and survival status in high- and low-risk groups, and the results demonstrated that the high-risk group presented a higher risk score and rates of death than the low-risk group (Figure 9D). The heatmap showed that the five gene had a higher expression level in the high-risk group than the low-risk group (Figure 9D). Similar results were found in the ICGC cohort (Figure 9E).

The HCC patients were split into high-risk and low-risk groups by the median risk score. In the training set, K–M analysis indicated that the patients with high-risk scores had worse overall survival than patients with low-risk scores (HR = 2.52, 95% CI:1.78-3.57, log-rank test p < 0.001) (Figure 10A). In the ICGC dataset, the overall survival of patients with low-risk scores was significantly longer than that of patients with high-risk scores (HR = 7.45, 95% CI: 4.09-13.59, log-rank test p < 0.001) (Figure 10B). To further understand the predictive ability at 0.5, 1, 2, 3, and 5 years of the signature, we used time–ROC analysis in TCGA (Figure 10C) and ICGC (Figure 10D) cohorts. The AUC of the signature in TCGA cohort was 0.764, 0.768, 0.712, 0.722, and 0.702 at 0.5, 1, 2, 3, and 5 years, respectively (Figure 9C). The AUC in the ICGC cohort was 0.725, 0.746, 0.801, and 0.823, 0.785 at 0.5, 1, 2, 3, and 5 years (Figure 10D). These suggest that the risk model provided a reliable predictive efficiency in HCC prognosis. Univariate Cox analysis was performed on all factors, and multivariate Cox analysis was performed with variates which had p < 0.05 of the univariate Cox analyses in TCGA (Figure 10E) and ICGC (Figure 10F) cohorts. The results confirmed that the clinical stage and risk score were stronger independent poor prognostic factors in both cohorts. To further determine the robustness of the risk mode, we performed an additional analysis to predict the overall survival of HCC patients in an independent GEO cohort (GSE14520). Similarly, K–M curves showed that the patient with high-risk had worse overall survival (Supplementary Figure S1A). The AUC of the signature was 0.73, 0.671, 0.610, 0.581, and 0.606 at 0.5, 1, 2, 3, and 5 years, respectively (Supplementary Figure S1B). The univariate and multivariate Cox regression analyses demonstrated that the risk score was an independent risk factor for the overall survival of HCC (Supplementary Figure S1C). This result was validated by another independent cohort CPTAC in the protein level, which completely echoed the result that high-risk patients appeared to present a poorer overall survival than those with low-risk scores. The AUC of the signature in the CPTAC cohort is 0.606 at 0.5 year, 0.706 at 1 year, 0.662 at 2 years, 0.672 at 3 years, and 0.702 at 5 years. Univariate and multivariate Cox analyses also demonstrated that protein risk score was an independent prognostic indicator in the protein level.

Considering the complexity of the relationship between ferroptosis and HCC, we illustrated the workflow of risk score construction with the alluvial diagram (Figure 11A). These results indicated that cluster A was linked to a lower risk score, whereas B and C clusters exhibited higher risk scores (Figure 11A). We also explored the relationship between the risk score and HCC subtype and found that the A subtype was associated with the lowest risk score, and C subtype was associated with the highest risk score (Figure 11B). The relationship between the risk score and alive or dead status showed that the patients with dead status were associated with a higher risk score than those with alive status (Figure 11C). The relationship between the risk score and difference grades showed that the patients with a higher grade appeared to be associated with a higher risk score (Figure 11D). The relationship between the risk score and difference in tumor stages showed that the patients with higher difference in tumor stages appeared to be associated with a higher risk score (Figure 11E). The relationship between the risk score and difference TNM stages showed that the patients with higher difference in TNM stage appeared to be associated with a higher risk score (Figure 11F). The alive patients were confirmed with low-risk score than those with high-risk score (alive vs dead: 77 vs 23% in the low-risk score group and 56 vs 44% in the high-risk score group; Figure 11G). The lower grade was confirmed in patients with low-risk score compared to those with high-risk score (G1 vs G2 vs G3 vs G4: 16 vs 55% vs 27 vs 2% in the low-risk score group, 11 vs 43% vs 41 vs 5% in the high-risk score group; Figure 11H). A consistent result was also observed in the tumor stage (Figure 11I) and TNM stage (Figure 11J). Collectively, these results revealed that the risk score is related to the tumor stage and correlates with the prognosis of HCC patients.

We built a comprehensive nomogram that could be applied in clinical practice based on the risk score and clinical stage (Figure 12A). Then, the calibration curves were performed to assess the prediction value of this nomogram (Figure 12B). The results demonstrated that the calibration curve for stage and risk score agreed well with the ideal observations. Furthermore, ROC analysis revealed that the combined nomogram showed the best performances in predicting the overall survival of HCC patients, compared to the risk score and clinical stage (AUC = 0.77 at 1 year, 0.76 at 3 years, and 0.74 at 5 years, Figure 12C). In addition, decision curve analysis indicated that the combined model achieved greater net benefit than the stage or risk score alone (Figure 12D).