The Prognostic Model Based on Tumor Cell Evolution Trajectory Reveals a Different Risk Group of Hepatocellular Carcinoma

Abstract

Hepatocellular carcinoma (HCC) is one of the leading causes of cancer-related death worldwide, and heterogeneity of HCC is the major barrier in improving patient outcome. To stratify HCC patients with different degrees of malignancy and provide precise treatment strategies, we reconstructed the tumor evolution trajectory with the help of scRNA-seq data and established a 30-gene prognostic model to identify the malignant state in HCC. Patients were divided into high-risk and low-risk groups. C-index and receiver operating characteristic (ROC) curve confirmed the excellent predictive value of this model. Downstream analysis revealed the underlying molecular and functional characteristics of this model, including significantly higher genomic instability and stronger proliferation/progression potential in the high-risk group. In summary, we established a novel prognostic model to overcome the barriers caused by HCC heterogeneity and provide the possibility of better clinical management for HCC patients to improve their survival outcomes.

Keywords: cell state transition; copy number aberration; genomic diversity; hepatocellular carcinoma; prognosis; single-cell transcriptomics; tumor evolution; tumor heterogeneity.

FIGURE 1

The whole process of data analysis.

FIGURE 2

Single-cell RNA sequencing (scRNA-seq) profiling of different malignant cell clusters. (A,B) The Uniform Manifold Approximation and Projection (UMAP) plot showing the annotation and color codes for cell types in hepatocellular carcinoma (HCC) (A). Cells were further shown in different color by patient origin (B). (C) The UMAP plot, showing only malignant cell clusters by Louvain algorithm. (D) With CopyKAT, malignant cell clusters were delineated into two subclones by single-cell copy number profiles inferred from scRNA-seq data. (E) Clustered heat maps of single HCC malignant cell copy number profiles in two major subclones.

FIGURE 3

Cells were sorted by progression from lower malignant state to higher malignant state. (A,B) Cell Trajectory performing the route of low- to high-malignant cells, which can serve as a model to describe malignant cell differences. (C) Expression levels for differentially expressed genes (rows), with cells (columns) shown in pseudo-time order. (D,E) ALB and Twist genes confirming the trusty of the cell progression trajectory model. (F,G) The Kaplan–Meier (K-M) analysis of the top 10 downregulated genes (root genes) and top 10 upregulated genes (end genes) from 1,000 differential genes group capturing the overall survival (OS) differences between low- and high-malignant cell groups. (H) Gene set variation analysis (GSVA) heat map showing the mainly differential signaling pathways between low- and high-malignant cell groups.

FIGURE 4

Establishing the 30-gene prognostic model with LASSO regression analysis. (A) LASSO regression analysis performed the frequency of different gene combination models and finally determined the 30-gene signature for OS prediction. (B) C-index of 30-gene prognostic model was 0.79 in TCGA training cohort, while 0.73 in ICGC validation cohort. (C) LASSO coefficient profiles of the gene features. (D) Ten-time cross-validation for tuning parameter selection in the LASSO model. (E,F) The risk score distribution and survival status distribution of 30-gene prognostic model in TCGA training cohort and ICGC validation cohort, and the heat map of gene expression are shown below with color, red (high) and green (low).

FIGURE 5

Prognostic performance of 30-gene signature in TCGA Training Cohort and ICGC Validation Cohort. (A,B) K–M survival curve for risk score in TCGA training cohort (A) and ICGC validation cohort (B). (C,D) Receiver operating characteristic (ROC) curve of the 30-gene prognostic model in TCGA cohort (C) and ICGC cohort (D). (E) Multivariate Cox regression analysis of clinical parameters and prognostic model for OS.

FIGURE 6

The analysis of genomic aberrations in high-risk group and low-risk group. (A) Recurrent copy number aberrations of high-risk group and low-risk group in TCGA cohort. Regions of recurrent copy number amplifications (red) and deletions (blue) were above and below baseline (0.0), respectively, in the targeted array were identified by GISTIC 2.0. (Red line represented GISTIC score of 0.3). (B,C) Oncoplot displaying the somatic landscape of high-risk group (B) and low-risk group (C). Genes were arranged according to their mutation frequency. The Y-axis was the gene name and the abscissa was the sample name. Different colors represented different mutation types. (D) Forest plot showed differentially mutated genes between high-risk group and low-risk group. The adjacent table included the number of samples in high-risk group and low-risk group with the mutations in the highlighted gene. The p-value indicated significance threshold: ***p < 0.001; **p < 0.01; Fisher’s exact test. (E) Co-bar plots showed the most recurrently mutated genes in high-risk group and low-risk group. (F) The distribution plot shows tumor mutation burden (TMB) distribution of different cancer types. Liver hepatocellular carcinoma (LIHC) patients were divided into low-risk group and high-risk group.

FIGURE 7

GSEA enrichment analysis. (A) The enrichment plot of upregulated gene sets in low-risk group. (B) The enrichment plot of downregulated gene sets in low-risk group.

FIGURE 8

Varied malignant cell subgroups contribute to the inter-tumor and intra-tumor heterogeneity of HCC. Hepatocyte-like tumor cells could progress to high plasticity tumor cells, accompanied by the inactivation of tumor suppressor pathways such as TP53, the disappearance of the inherent characteristics of hepatocytes, the enhancement of proliferation, invasion and metastasis ability, and the appearance of immune suppression.