Prognostic iron-metabolism signature robustly stratifies single-cell characteristics of hepatocellular carcinoma

Abstract

Cancer immunotherapy has shown to be a promising method in treating hepatocellular carcinoma (HCC), but suboptimal responses in patients are attributed to cellular and molecular heterogeneity. Iron metabolism-related genes (IRGs) are important in maintaining immune system homeostasis and have the potential to help develop new strategies for HCC treatment. Herein, we constructed and validated the iron-metabolism gene prognostic index (IPX) using univariate Cox proportional hazards regression and LASSO Cox regression analysis, successfully categorizing HCC patients into two groups with distinct survival risks. Then, we performed single-sample gene set enrichment analysis, weighted correlation network analysis, gene ontology enrichment analysis, cellular lineage analysis, and SCENIC analysis to reveal the key determinants underlying the ability of this model based on bulk and single-cell transcriptomic data. We identified several driver transcription factors specifically activated in specific malignant cell sub-populations to contribute to the adverse survival outcomes in the IPX-high subgroup. Within the tumor microenvironment (TME), T cells displayed significant diversity in their cellular characteristics and experienced changes in their developmental paths within distinct clusters identified by IPX. Interestingly, the proportion of Treg cells was increased in the high-risk group compared with the low-risk group. These results suggest that iron-metabolism could be involved in reshaping the TME, thereby disrupting the cell cycle of immune cells. This study utilized IRGs to construct a novel and reliable model, which can be used to assess the prognosis of patients with HCC and further clarify the molecular mechanisms of IRGs in HCC at single-cell resolution.

Keywords: Hepatocellular carcinoma (HCC); Immune infiltration; Iron metabolism-related gene prognostic index (IPX); Iron metabolism-related genes (IRGs); Personalized medicine.

Graphical abstract

Fig. 1

The schematic workflow of the study.

Fig. 2

Prognostic DIRGs could classify patients with HCC into three clusters with distinct prognostic outcomes. (A) Volcano plot showing DIRGs based on TCGA training cohort. (B) Heatmap showing similarity matrix of patients in TCGA training cohort derived from unsupervised consensus clustering. (C) Kaplan-Meier plot showing survival difference within three clusters. (D) Forest plot showing prognostic DIRGs associated with OS.

Fig. 3

Performance of the predictive score model (A-C) and the nomogram (D-F). (A, D) Kaplan-Meier survival curves for comparison of the overall survival rates between patients in the low-risk group and the high-risk group for the training cohort (upper), The 1,3,5-year ROC curve of predictive score model for the training cohort (lower, all P-value < 0.001). (B, E) Comparison of the overall survival rates between patients in the low-risk group and the high-risk group for internal testing cohort (upper), The 1,3,5-year ROC curve of predictive score model for internal testing cohort (lower, all P-value < 0.001). (C, F) Comparison of the overall survival rates between patients in the low-risk group and the high-risk group for external testing cohort (upper), The 1,3,5-year ROC curve of predictive score model the external testing cohort (lower, all P-value < 0.001).

Fig. 4

mRNA and protein expression levels of individual genes in DIRGs in different datasets. The expression differences of 6 of the 13 DIRGs were statistically significant, and the expression trends were consistent in the TCGA dataset (A-F), ICGC dataset (G-L), The Human Protein Atlas (M-Q), and 12 HCC clinical samples (using RT-qPCR) (R-W).

Fig. 5

Detailed characterization of malignant cells based on single-cell RNA sequencing data. (A) UMAP visualization of signature genes that were used to annotat major cell types. (B) UMAP visualization of integrated 6 major cell types. (C) UMAP visualization of hepatocytes that were clustered into 15 sub-clusters. (D) Stacked barplot showing the proportion of hepatocytes cell types in each site. (E) Violin plot showing the risk scores of hepatocytes across different HCC stages. (F) Violin plot showing the risk scores of hepatocytes across different sites. (G) Violin plot showing the risk scores of hepatocytes across different hepatocyte sub-clusters. (H) Module–trait relationships in all hepatocyte sub-clusters. (I) Relatedness network of genes in the red module. (J) The red module membership and gene significance for the C4 sub-cluster. (K) The red module membership and gene significance for the C5 sub-cluster. (L) Heatmap showing the regulon modules using regulon connection specificity index (CSI) matrix with representative regulators. (M) Rank plot showing the transcription factors of C4 (Left) and C5 (Right) sub-cluster based on regulon specificity score (RSS).

Fig. 6

Discriminating T-cell populations based on single-cell RNA sequencing data. (A) Stacked barplot showing the proportion of CD4⁺ T cells between the high-score and low-score populations. (B) The trajectory of CD8⁺ T cells inferred by Monocle2. (C) Pseudotime analysis showed six CD8⁺ T subpopulations on each branch. (D) Pseudotime analysis showed high-score and low-score CD8⁺ T cells on each branch. (E) Stacked barplot showing the proportion of CD8⁺ T cells between the high-score and low-score populations. (F) The trajectory of CD4⁺ T cells inferred by Monocle2. (G) Pseudotime analysis showed six CD4⁺ T subpopulations on each branch. (H) Pseudotime analysis showed high-score and low-score CD4⁺ T cells on each branch. (I) Pseudotime plot showing distinct developmental trajectories of CD4⁺ T cell and CD8⁺ T cell for high-score and low-score populations. (J) Heatmap for the dynamic top genes along pseudotime for high-score and low-score populations for CD4⁺ and CD8⁺ T cells.

Fig. 7

Enrichment analysis showed distinct biological processes between high-score and low-score clusters. (A) Biological pathways that are associated with significantly branch-dependent genes for high-score. (B) Biological pathways that are associated with significantly branch-dependent genes for low-score. (C) Single sample gene set enrichment analysis results for comparing the expression of the dynamic genes during branch evolution for the “Cell cycle” term, the “Positive regulation of cell activation” term, the “Treg” term, and the “T exhaustion” term. * P-value < 0.05, ** P-value < 0.01, *** P-value < 0.001 and **** P-value < 0.0001.