The normalized gene-level RNA-Seq data and clinical information for 364 patients LIHC-TCGA cohorts were downloaded from UCSC Xena¹ with R package UCSC Xena Tools (Wang and Liu, 2019). To obtain 258 patients LIRI-JP validation set, RNA-seq data, and related clinic pathological data were downloaded from the ICGC website² (Zhang J. et al., 2019). The scRNA-seq barcode sequences and raw gene expression matrix were downloaded from the CNP0000650 (Sun et al., 2021). Mutation data that contained somatic variants were stored in Mutation Annotation Format (MAF) form and were downloaded from Genomic Data Commons (GDC).³

To verify the trend of this tumor evolution trajectory, Kaplan--Meier (K--M) analysis was performed. The top 10 end-genes were extracted and the potential prognostic significance of these genes was assessed with the LIHC data from GEPIA2.⁴ The K–M survival curves were also generated to graphically demonstrate the OS to the high-risk group and low-risk group, which were stratified by the tumor evolution signature. The R package called “survminer” was utilized to perform the survival analysis, and the optimal cutoff was ascertained by the “surv_cutpoint” function.

Mutation comment file (MAF) of TCGA-LIHC cohort was downloaded from the GDC client. Differential analysis and visualization of somatic mutations was done using Maftools (The Cancer Genome Atlas Research Network et al., 2013; Mayakonda et al., 2018). This difference between high- and low-risk group was detected using function “mafComapre,” which performs Fisher’s exact test on all genes between two groups to detect differentially mutated genes.

Composite copy number profiles were generated to highlight differences between high- and low- risk group. Segment file of TCGA-LIHC cohort was downloaded from FIREHOSE and samples were further divided into high- and low-risk groups. Then we ran the GISTIC 2.0 pipeline to generate discrete copy number data file. Chromosomes reference objects were from the “BSgenome.Hsapiens.UCSC.hg19” R package.

As in the previous study, a non-synonymous mutation from the TCGA database was used as the raw mutation count, and it was divided by 38 MB to quantify TMB (Chalmers et al., 2017). The samples were sorted according to the value of the median TMB from low to high.

Gene Enrichment Analysis (GSEA) was further used to investigate the functional enrichment of risk score-associated genes using the R package “clusterProfiler” (Yu et al., 2012). The Benjamini–Hochberg method was used to adjust nominal p-values (false discovery rate, FDR) for multiple testing.

The Maftools package was used to illustrate the respective mutation profiling of the two risk group levels by waterfall plot, and differentially mutated genes were identified by using the “mafCompare” function where genes mutated in greater than 5% of LIHC samples in the TCGA cohort were considered (Mayakonda et al., 2018).

Student’s t-test was conducted to make statistical comparison. The “pheatmap” R package was applied to generate heatmaps. Survival analysis was completed using Kaplan–Meier method, and the prediction performance of the risk model was evaluated using receiver operating characteristic (ROC) via “time-ROC” R package. Multivariate COX regression analyses were used to investigate the prognostic value of risk-score. Hazard ratio (HR) and 95% confidence intervals (CI) for each variable were also calculated where needed. A value of p < 0.05 was defined as statistically significant difference. All of our analyses were conducted using R software version 4.0.2.⁵

The whole process of data analysis is depicted in Figure 1.

Unsupervised dimensionality reduction and graph-based clustering analysis were performed with the data from CNP0000650, and 24 clusters (Figure 2A) were visualized by the UMAP method (Becht et al., 2019; Sun et al., 2021). The immune cells mainly consisted of myeloid-derived cells, T cells, B cells, plasma cells, and natural killer (NK) cells, while non-immune cells included endothelial cells, hepatic stellate cells, apparently normal epithelial cells, and HCC malignant cells. To contrast the difference among different patients, we classified the cells by patient origin (Figure 2B), and the result showed that tumor cells contained obvious heterogeneity, while non-tumor cells kept homogeneous, proving that the differences between tumor cell clusters are mainly due to the tumor heterogeneity, rather than batch effects between samples. No normal liver cells were detected, likely because of the technical limitations, resulting in no comparison between normal and malignant liver cells.

We further extracted the varied tumor cell clusters for the analysis of tumor heterogeneity (Figure 2C), and two major sub-clones were defined by the clustered heat maps of single cell copy number profiles (Figure 2D). Compared with sub-clone 2 (green), the heatmap showed that sub-clone 1 (red) contained more CNAs, implying that sub-clone 1 might be more malignant (Figure 2E).

To determine the relationship between malignant cell clusters, we performed single-cell trajectory analysis with scRNA-seq data using Monocle (Qiu et al., 2017). As is well known, genomic mutations accumulate over time in the process of tumor evolution (Turajlic et al., 2019; Craig et al., 2020; Losic et al., 2020). Thus, we defined the cell cluster with a lower CNA burden as the root, while the cell cluster with a higher CNA burden was defined as the end of the trajectory. We noticed that sub-clone 2, comprising cells with obvious liver characteristics, was concentrated at the beginning of the trajectory, while sub-clone 1, comprising cells with less specificity of origins, was concentrated at the end of the trajectory, indicating that the trajectory model fits the process of tumor evolution well (Figures 3A,B).

During the process of transition of tumor cells from a lower to higher malignant state, some genes are silenced, while others become newly active. We used Tradeseq, a powerful generalized additive model framework based on the negative binomial distribution, to interpret the within lineage differential expression (Figure 3C and Supplementary Table 1). To verify the malignant trend of this trajectory, we extracted the top 10 downregulated genes (root genes) and top 10 upregulated genes (end genes) for Kaplan–Meier (K-M) analysis. The high expression of the top 10 root genes represented a benign prognosis, while the high expression of the top 10 end genes represented a poor prognosis (Figures 3F,G). Furthermore, along the trajectory, liver characteristics such as ALB (Figure 3D) were gradually lost, while stemness and malignant marker genes such as TWIST1 (Figure 3E), gradually increased, suggesting that the malignant state was advancing.

The gene set variation analysis (GSVA; Hänzelmann et al., 2013) was used to further analyze the underlying biological processes along the trajectory. In the Molecular Signature Database (MSigDB) “hallmark” collection of major biological categories (Liberzon et al., 2015), the upregulated genes of sub-clone 1 were enriched in the tumor-promoting pathway (“Wnt/β-catenin signaling”) and proliferation pathway (“G2M checkpoint,” “E2F Targets,” “MYC Targets”), while the downregulated genes of sub-clone 1 were enriched in the tumor-suppressor pathway (“P53 pathway”) and essential liver function pathway (“Complement,” “Fatty acid metabolism,” “Adipogenesis”) (Figure 3H), which were consistent with the characteristics of cells that progressed from the lower malignant state to the higher malignant state (Maley et al., 2017; Losic et al., 2020; Marjanovic et al., 2020; Llovet et al., 2021).

Although the top 10 genes had a certain predictive effect, we preferred optimizing gene combination to obtain a better prognostic model. With the selection criteria of p < 0.01, the intersection of univariate Cox regression analysis and K-M analysis identified 200 credibly survival-related genes. We used TCGA data as the training cohort and ICGC data as the external validation cohort (The Cancer Genome Atlas Research Network et al., 2013; Zhang J. et al., 2019). Lasso-penalized Cox analysis was subsequently performed 1,000 times in the TCGA training cohort with 10-fold cross-validation to evaluate and eliminate variables that contributed less to the model, and a 30-gene signature with the most powerful predictive features were selected (Figures 4A,C,D). To validate the credibility of this model, C-index was assessed in the TCGA training cohort and ICGC validation cohort, which was confirmed as 0.79 and 0.73, respectively (Figure 4B), suggesting that our model had favorable efficacy for predicting prognosis (Harrell, 1982). Based on the 30-gene prognostic model, TCGA and ICGC samples were clustered into high-risk and low-risk groups, and the OS time of patients in the high-risk group was remarkably decreased (Figures 4E,F).

The K–M analysis and time-dependent ROC was used to assess the prognostic capacity of the 30-gene prognostic model in the TCGA cohort and ICGC cohort, respectively. The K–M analysis illustrated that patients in the low-risk group had significantly longer OS than those in the high-risk group, both in the TCGA cohort (Figure 5A) and the ICGC cohort (Figure 5B). The area under the ROC curve (AUC) for the 1-, 3-, and 5-year OS was 0.843, 0.848, and 0.824 in the TCGA cohort, while it was 0.77, 0.796, and 0.774 in the ICGC cohort (Figures 5C,D), indicating that this 30-gene prognostic model had high sensitivity and specificity for survival prediction.

Gender, age, stage, vascular invasion, bile duct invasion, fibrosis, and the risk score of the prognostic model were included in the multivariate Cox regression model, and the risk score was revealed to be independent predictor for OS, showing splendid predictive performance ability, with HR: 6.432, 95% CI: 4.095–10.103, p < 0.001 in TCGA cohort, and HR: 3.751, 95% CI: 1.419–9.918, p < 0.001 in ICGC cohort. Taken together, the 30-gene prognostic model was completely reliable for the precise prediction of OS in HCC (Figure 5E).

Increasing mutation frequency is a typical feature of human cancer. To identify the divergence of genomic aberrations between the high-risk and low-risk groups in the TCGA database, CNAs data were downloaded from the GDC portal and analyzed with GISTIC 2.0, with which the high-risk group had an obviously higher genomic aberration burden than the low-risk group (Figure 6A). Meanwhile, the top 20 genes with high genomic mutation frequency in the high-risk and low-risk groups were constructed by Maftools (Mayakonda et al., 2018; Figures 6B,C). To analyze the discrepancy between the high-risk and low-risk groups, the differentially mutated gene type and frequency were compared by Fisher’s exact tests. The results showed three significantly differential genes—TP53 (47 versus 19%), OBSCN (18 versus 5%), and RB1 (11 versus 2%) (Figure 6D). The six genes most recurrently mutated were TP53 (47 versus 19%), TTN (31 versus 27%), CTNNB1 (26 versus 28%), MUC16 (19 versus 16%), OBSCN (18 versus 5%), and ALB (12 versus 13%) (Figure 6E and Supplementary Table 2). These findings suggested that some mutated genes such as TP53 and OBSCN that were notably different compared with the high-risk and low-risk groups could be related to the malignant progression and can continuously accumulate mutations over time; whereas, the others such as CTNNB1 and ALB, which remained stable from the low-risk state to high-risk state, likely contribute to the essential neoplastic process rather than malignant progression.

The tumor mutation burden (TMB) is also considered an essential factor impacting on the occurrence and progression of the tumor. The distribution plot shows TMB distribution of different cancer types (Figure 6F). The three types including low-risk HCC, all HCC samples (liver hepatocellular carcinoma, LIHC), and high-risk HCC were apparently distinct from each other, and the TMB gradually increased from low-risk type to high-risk type, which suggested that our model had excellent distinguishing capability.

To explore the underlying molecular mechanisms of this prognostic model, we conducted GSEA to compare the low-risk group with the high-risk group in TCGA cohorts (Yu et al., 2012). In the MsigDB “hallmark” collection of major biological categories, proliferative signaling pathways (“E2F targets,” “G2M checkpoint,” “KARS signaling”), and the invasion and metastasis-related signaling pathways (“EMT” and “myogenesis”) were dramatically increased in the high-risk group (Figure 7B). This was consistent with the data at the single-cell level. Notably, the low-risk group was enriched in inflammation-related gene pathways such as “inflammatory response,” “interferon-α response,” and “interferon-γ response,” which suggested that the secretion of inflammatory factors might originate from the tumor cells, indicating a strong proinflammation potential in the early tumorigenesis stage (Figure 7A). A similar phenomenon had been reported in melanoma, breast cancer, and colorectal cancer, where the expression of immune-related genes was also presented by tumor cells and could likely be an independent influence on their prognostic differences (Sconocchia et al., 2014; Forero et al., 2016; Buetow et al., 2019; McCaw et al., 2019; Jin et al., 2020).