Identification of two gene expression-based clusters in HCC
The discovery cohort consisted of 377 HCC profiles profiled by TCGA, and the validation cohort consisted of 247 HCC profiles profiled by Affymetrix (Supplementary Table 1). To identify the candidate genes, we developed a strategy. Firstly, genes with a mean expression < 1 in the tumor tissue were excluded. Next, the genes were ranked by their coefficient of variation (Supplementary Table 2). The top 10% of most variable genes were used as input for hierarchical clustering (HC). The discovery cohort was divided into two clusters (Fig. 1A), with the larger of which (n = 209) was designated category A (CA) and the smaller cluster (n = 162) was designated category B (CB) (Supplementary Table 3). In order to confirm the robustness of this strategy we randomly selected 70% of the samples as the training set, and used the aforementioned top 10% genes with the highest variation and SVM algorithm to train the model. Then, based on this model, we divided the validation set. Finally, we compared it with the original labels and plotted the ROC curve (AUC = 0.869Supplementary Fig. 1). Patients with CA tumors showed a higher prevalence among females, younger age, higher AJCC stage, and worse overall survival (OS) (Hazard ratio (HR) = 1.500, log-rank test p = 0.0266, 95% CI 1.048–2.148; Fig. 1B, Supplementary Table 4). The multivariate Cox regression analysis of clinical factors revealed that CA and AJCC stage III/IV were strong prognostic factor for patient survival, while AJCC stage II was not (p = 0.041, p < 0.001, p = 0.014, respectively; Fig. 1C).
Identification of two gene expression-based clusters in HCC. A, HC analysis identified two molecular clusters at the transcriptome level in TCGA-LIHC. B, Kaplan-Meier curves compared OS between CA and CB in TCGA-LIHC. C, Forest plot showed the results of multivariate Cox regression analysis of molecular subtyping, AJCC staging system, ethnicity, and age.
Comparison of AJCC staging system and molecular subtyping
To investigate whether the molecular subtyping could be comparable to AJCC staging system in terms of prognostic prediction, we conducted a series of survival analysis. Consistent with the previous results, patients with stage III/IV exhibited outcomes than those with stage I (HR = 2.676, log-rank test p < 0.001, 95% CI 1.754–4.083; HR = 5.496, log-rank test p = 0.00452, 95% CI 1.695–17.821; Fig. 2A). However, there was no significant difference in survival between stage II and stage I. We hypothesized that the lack of significance might be attributed to the different versions of AJCC staging system. To standardize the criteria, we conducted a restaging of cases. TCGA provides a detailed record of the TNM staging for HCC, as well as the AJCC version. According to the above information, we translated the different versions of AJCC into the 8th edition. After restaging, we did not observe more pronounced survival differences, yet (Fig. 2B). This situation may be due to the lack of “regional lymph node” and “metastasis” information, which caused incorrect judgement (Supplementary Table 5). In view of the situation above, the samples with ambiguous judgement of stage were excluded, resulting in the acquisition of a more discriminating curve (Fig. 2C). These findings suggested that AJCC 8th staging system exhibited superior prognosis capability compared to previous versions. Despite its strong predictive ability, there was still no statistically significant differentiation in prognosis between stage I and stage II patients. While, our study revealed that molecular subtyping showed more promising performance in early-stage patients (HR = 1.785, log-rank test p = 0.0818, 95% CI 0.930–3.429; Fig. 2D). As a result, we devised a new stratification strategy, reclassifying patients based on both AJCC 8th stage and molecular subtyping. The patients were reclassified into different groups according to AJCC 8th stage and molecular subtyping. This approach unveiled significantly worse outcomes for patients in stage II and category A subtype compared to stage I and category B subtype (HR = 2.542, log-rank test p = 0.0297, 95% CI 1.096–5.894; Fig. 2E). Notably, all stage IV patients belonged to category A, prompting further exploration of the applicability of molecular subtyping in stage III (Supplementary Table 4). Although a trend was observed with shorter overall survival (OS) in category A compared to category B patients, it did not achieve statistical significance (Supplementary Fig. 2). Further analysis revealed that patients in stage III/IV with category A experienced worse prognosis than others (HR = 3.362, log-rank test p < 0.001, 95% CI 2.066–5.472; Fig. 2F).
Comparison of AJCC staging system and molecular subtyping. A, Kaplan-Meier curves compared OS among patients with different AJCC stage. B, Kaplan-Meier curves compared OS among patients with adjusted AJCC 8th stage. C, Kaplan-Meier curves compared OS among patients with confirmed AJCC 8th stage. D, Kaplan-Meier curves compared OS between CA and CB in early stage. E, Kaplan-Meier curves compared OS between patients with different groups defined by stage and molecular subtyping in early stage. F, Kaplan-Meier curves compared OS between patients with different groups defined by stage and molecular subtyping in all stage.
Cluster allocation of validation cohorts
To assess the predictive efficacy of our molecular subtyping, we utilized the same candidate genes for HC on the GSE14520 dataset, resulting in the categorization of patients into CA (n = 145) and CB (n = 102) groups (Supplementary Fig. 3, Supplementary Table 6). Notably, individuals with CA tumors exhibited a younger age, higher CLIP stage, elevated AFP levels, and poorer clinical outcomes. We observed significant overlap between our CA cluster and GSE14520’s “high-risk metastasis signature”, as well as between the CB cluster and GSE14520’s “low-risk metastasis signature” (p < 0.0001) (Supplementary Table 7). As expected, AJCC staging system demonstrated outstanding predictive performance in terms of OS and recurrent-free survival (RFS) (Supplementary Fig. 4–5). Similarly, patients with CA tumors experienced inferior OS (HR = 1.672, log-rank test p = 0.0168, 95% CI 1.097–2.549; Fig. 3A) and RFS (HR = 1.539, log-rank test p = 0.0152, 95% CI 1.087–2.178; Fig. 3B) compared to patients with CB tumors. However, our molecular subtyping on its own was insufficient to independently predict the prognosis of patients in the validation cohorts (Fig. 3C-D). Notably, the combined use of molecular subtyping and AJCC staging system outperformed using AJCC staging system alone (Fig. 3E-F). Even in stage III, patients with CA tumors exhibited worse outcomes than those with CB tumors (Supplementary Fig. 6–7). Furthermore, it was observed that patients with stage III and CA tumors had the worst prognosis (Fig. 3G-H).
Cluster allocation of validation cohorts. A-B, Kaplan-Meier curves compared OS (A) and RFS (B) between CA and CB in all stage. C-D, Kaplan-Meier curves compared OS (C) and RFS (D) between CA and CB in early stage. E-F, Kaplan-Meier curves compared OS (E) and RFS (F) between patients with different groups defined by stage and molecular subtyping in early stage. G-H, Kaplan-Meier curves compared OS (G) and RFS (H) between patients with different groups defined by stage and molecular subtyping in all stage.
Multi-omics analyses of prognosis-related clusters
To explore potential differences between CA and CB tumors at the multi-omics level, we obtained simple nucleotide variation (SNV) (n = 362), copy number variation (CNV) (n = 358), and DNA methylation (DNAm) (n = 371) data that overlapped with transcriptome profiling data for LIHC from TCGA.
Firstly, we evaluated the general landscape of SNV and CNV at the genomic level (Fig. 4A, Supplementary Fig. 8–9). Notably, 75 genes (defined as DMUTG) exhibited significant differences in mutation frequency between CA and CB at p < 0.05 (Supplementary Table 8). CA tumors demonstrated a prevalence of mutations in TP53, CHD4, MLXIPL, IRF2, MTMR4, BAP1, and TLR8, which are associated with molecular functions such as DNA transcription factor binding, chromatin binding, and protein phosphatase binding. In contrast, CB tumors showed a predominance of mutations in CTNNB1, ATP10A, ATP2B4, RASA3, ACVR2A, ALOX12B, ITGA9, NUP153, and others, associated with molecular functions like metal ion binding, ATP binding, catenin binding, and protein kinase activity. Notably, TP53 mutations were mutually exclusive with CTNNB1 mutations significantly in HCC (p = 0.0357), indicating different mutation preferences for CA and CB. Furthermore, we computed mutational signatures to explore the connection between molecular subtyping and specific mutational processes based on Catalogue Of Somatic Mutation In Cancer (COSMIC) (Supplementary Fig. 10). Our analysis revealed that signature 16, signature 4, signature 3, signature 24, signature 1, and signature 22, which were commonly identified in liver cancers, accounted for half of mutational signatures. Interestingly, CA showed a predominance of signature 1, associated with an endogenous mutational process initiated by spontaneous deamination of 5-methylcytosine, while CB tumors showed a prevalence of signature 22, associated with exposures to aristolochic acid (Fig. 4B).
Based on segmented copy number data from TCGA, we identified 67 focal candidate copy number alterations (CNAs), involving 1191 genes at q < 0.1 (Supplementary Table 9). Specifically, when considering CA and CB separately, we found more CNAs and genes, with 46 CNAs and 1679 genes for CA and 48 CNAs and 1839 genes at q < 0.1 (Supplementary Table 10–11, Supplementary Fig. 11–14). These results suggested that mixture of CA and CB may obscure the true signals of different subgroups. The genes amplified in CA and deleted in CB or amplified in CB and deleted in CA were defined as differential copy number alterations genes (defined as DCNAG, Fig. 4C). At the genome level, we identified 75 DMUTG and 61 DCNAG (defined as DGG together) between CA and CB.
Furthermore, CA and CB exhibited significant differences in the performance of transcriptome and epigenome. We identified 3542 differentially expressed genes (defined as DEG) and 2101 differentially methylated loci at fdr < 0.05 and |log2 fold change| > 1. These loci involved 1029 genes. Among them, 1026 genes had the same differentially methylated state (different probes targeting the same gene all indicate that the gene is hypermethylation or hypomethylation, defined as DMG) between CA and CB (Supplementary Table 12–13, Fig. 4D-E). Summarizing these finding at multi-omics landscape, we discovered that PTPRT and ARHGEF4, which were overexpressed in CA, altered simultaneously. However, their specific contributions to the prognosis of HCC remained unknown.
Multi-omics analyses of prognosis-related clusters. A, Plot of TGCA-LIHC. The center plot represented mutation type and frequencies for 30 high-frequency genes. The histogram above the main plot represented tumour mutation burden (TMB, mutations/megabase) per sample. The histogram to the left displayed the statistical significance of mutation frequency between CA and CB. The histogram to the right represented the proportion of different mutation types. The heatmap below the main plot represented the clinical information of patients. B, Boxplot compared the proportion of different mutational signatures between CA and CB. C, Venn plot showed CNV between CA and CB. D-E, Volcano plot showed differential expressed genes (D) and differential methylated genes (E) between CA and CB. F, Venn plot showed differences between CA and CB at the multi-omics level.
Association between biological pathways, TIME and molecular subtyping
To explore the impact of transcriptome changes on biological functions, GSEA was utilized. We discovered 18 activated and 22 suppressed pathways in CA tumors, as evidenced by the large differences between CA and CB (Supplementary Table 14, Fig. 5A). Notably, CA tumors exhibited aberrantly excited cytokine receptor interaction and cell proliferation, while CB tumors were characterized by unusual activation of lipid metabolism, amino acid metabolism, and glycometabolism.
Given that the liver plays a crucial role in both immunity and immune tolerance, the immune microenvironment not only profoundly impacts disease progression and immunotherapy but also affects patient survival in HCC. In this research, we assessed the immune infiltration in HCC using the TCGA-LIHC dataset (Supplementary Fig. 15–16) and performed clustering analysis on the results. We found that our molecular subtyping differed significantly from the immune subtype predicted by CIBERSORT and ssGSEA (Fig. 5B-C). However, based on the evaluation results of ssGSEA, we found that CA tumors were highly infiltrated by immune cells (Supplementary Fig. 17). Interestingly, we found that 47 patients with CA tumors (203 in all) and 85 patients with CB tumors (159 in all) suffered from four gene mutations (CTNNB1, AXIN1, APC, and NCOR1). These four genes involved in the oncogenic Wnt pathway, which could alter immune infiltration in HCC 32.
Association between biological pathways, TIME and molecular subtyping. A, Dot plot of GSEA between CA and CB. B-C, ROC curves for the CIBERSORT (B) and ssGSEA (C) to predict molecular subtyping.