The research scheme is illustrated in Figure 2 . Further details on analytic processes and experimental methods are provided in the Supplementary Materials .

Bulk RNA sequencing (RNA-Seq) data were obtained from the Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC) and the International Cancer Genome Consortium (ICGC-LIRI-JP) projects. A total of 89 feature scores (FSs) were calculated, including gene set variation analysis (GSVA) scores (28), curated scores, and immune cell infiltration metrics. GSVA scores encompassed 50 hallmark gene sets and another 11 crucial gene sets ( Supplementary Table S1 ). Six curated scores—mRNAsi, TIDE, estimate score, immune score, stromal score, and tumor purity—were calculated based on previous studies (22, 25, 26). Immune cell infiltration was estimated using the R package “immunedeconv”, which involves six algorithms: Cibersort, Quantiseq, Timer, MCPcounter, XCell, and EPIC (21). Specifically, the proportions of 22 immune cell types calculated by the Cibersort algorithm were enrolled in the 89 FSs above.

Differential, prognostic, and Venn analyses were performed to identify consistently significant FSs in HCC. These FSs were then used for the unsupervised clustering analysis (29). The characteristics of clusters were analyzed, including key signal pathways, immune evaluation, pseudo-temporal trajectory, single nucleotide variation (SNV), drug sensitivity, and single-cell RNA sequencing (scRNA-Seq). Moreover, to make out the association among diverse classification criteria, we compared the proportion of different stages and MFP subtypes among clusters, respectively (30).

Cluster-specific marker genes were identified employing the “FindAllMarkers” function from the “Seurat” package (31). A classifier was developed using the random forest (RF) algorithm based on these marker genes in the TCGA-LIHC training cohort and externally validated with the ICGC-LIRI-JP cohort. To validate the reliability of our classifier, characteristics of the predicted classifications in the external cohort were compared with those of the clusters in the TCGA-LIHC cohort. Additionally, datasets GSE202069 and PRJEB34724 were used to validate the responses to anti-PD1 therapy based on our classifier (12, 32, 33).

A novel prognostic nomogram was developed using cluster-specific marker genes in the TCGA-LIHC training cohort and externally validated with the ICGC-LIRI-JP cohort.

The practical software was developed for individualized predictions of molecular classifications and prognoses using external transcriptomic data.

Significant therapeutic targets were identified through Venn analysis, differential expression analysis, and prognostic analysis of genes from several prognostic models. Lentiviral transduction was used to establish overexpressing Huh7 cell lines for the target gene and its negative control. The biological roles of the gene were evaluated by means of cell proliferation, colony formation, transwell migration, apoptosis, cell cycle assays, and tumorigenicity assays in nude mice.

All animal procedures were approved by the Ethics Committee of Zhejiang Cancer Hospital and conducted in accordance with the ethical guidelines.

After calculating the 89 FSs, differential and prognostic analyses were performed to identify the most critical FSs in the TCGA and the ICGC cohorts. The differential analysis results are presented in the heat maps ( Figures 3A, B ) and Supplementary Table S2 , while the Kaplan-Meier survival analysis results are provided in Supplementary Table S3 . The forest plots of the univariable Cox regression analyses are shown in Supplementary Figure S1 . The results of the Kaplan-Meier survival analysis were incorporated into the subsequent screening. Finally, Venn analysis for the differential and the prognostic FSs identified 37 critical FSs closely related to HCC biological behavior ( Figure 3C ).

Using the 37 critical FSs, we performed unsupervised clustering analysis on 374 tumor samples from the TCGA-LIHC cohort, identifying four clusters: C1, C2, C3, and C4 ( Figures 3D–F ).

Significant survival differences were observed across the clusters. Patients in C3 had the most favorable prognosis, with a median overall survival (OS) of 81.9 months, followed by C2 (median OS: 71 months), C1 (median OS: 45.7 months), and C4 (median OS: 30 months) (P < 0.001, Figure 4A ).

The results showed that C3 had the highest proportion of stage I tumors (58.14%) and the lowest proportion of stage III (12.79%), while C4 had the highest proportion of stage III tumors (35.44%) and the lowest proportion of stage I (31.65%). C1 was predominantly MFP-IE (37.50%), whereas C2 and C4 were mainly MFP-D (64.29% and 72.15%, respectively). C3 showed the highest frequency of MFP-F (47.67%) ( Figure 4B ).

To investigate the evolutionary relationships among clusters, we performed a pseudo-temporal trajectory analysis. As illustrated in Figure 4C , the lower-right subpopulation, representing normal samples, was designated as the starting point, while the top-left subpopulation marked the endpoint. The analysis revealed that cluster C3 was closest to normal tissue, followed by C2, C1, and C4, which aligns closely with their prognostic characteristics ( Figure 4D ).

Significant differences in gene expression and FSs were observed among clusters. The expression of MKI67 and AFP was significantly upregulated in C1 and C4 than that of C2 and C3. Compared to the normal samples, samples with upregulated AFP in HCC tumors mostly occurred in C1 and C4, not C2 and C3 ( Supplementary Figure S2 ; Supplementary Table S4 ). The significant feature scores (FSs) in the heat map can be mainly categorized into three parts: stemness proliferation, microenvironment (immune/stroma), and metabolism. In terms of stemness proliferation, notable FSs include the mRNAsi score, MYC targets, E2F targets, G2M checkpoint, PI3K AKT mTOR pathway, and DNA repair activity. The microenvironment category encompasses the estimate score, stromal score, immune checkpoint score, immune score, complement score, immunogenic cell death score, tumor purity, and TIDE score. The metabolism category includes various metabolic-related pathway scores. These categories help to delineate the roles and impacts of FSs within the context of tumor biology.

C1 was characterized by upregulated activities in MYC targets, E2F targets, G2M checkpoint, and PI3K AKT mTOR. In addition, C1 had a high estimate score, hot immunity (such as immune checkpoint score, immune score, complement score, and immunogenic cell death score), and high disulfidptosis, while its tumor purity and TIDE score were low, as well as hypometabolism; C2 was characterized by downregulated PI3K AKT mTOR activity, cold immunity, high tumor purity, and hypermetabolism; C3 exhibited the lowest stemness, high estimate score and stromal score, hot immunity, low tumor purity, low DNA repair activity, and hypermetabolism; C4 showed the highest stemness, high G2M, high DNA repair, but low PI3K AKT mTOR activity, cold immunity, low estimate score and stromal score, high tumor purity, and hypometabolism ( Figure 4E ; Supplementary Table S5 ).

GSEA revealed significant pathway alterations among the clusters. In C1, there was upregulation in pathways such as epithelial-mesenchymal transition (EMT), interferon response, and antigen processing, while pathways like adipogenesis, fatty acid metabolism, and oxidative phosphorylation were downregulated. Conversely, C2 exhibited the opposite trends in these pathways compared to C1. C3 showed upregulation in EMT, interferon response, coagulation, and complement pathways, with downregulation observed in the unfolded protein response, G2M checkpoint, E2F targets, and MYC targets. In C4, significant upregulation was observed in the G2M checkpoint, E2F targets, ribosome, and oxidative phosphorylation pathways, while its peroxisome, immune, and metabolic pathways were significantly downregulated ( Supplementary Figure S3 ; Supplementary Table S6 ).

We characterized the immune landscape among the four clusters to infer their responses to immunotherapy. Immune and stromal cell infiltration was evaluated using six computational methods: XCell, MCPcounter, Timer, Quantiseq, EPIC, and Cibersort ( Supplementary Figure S4A ). Overall, clusters C1 and C3 had higher levels of immune and stromal cells compared to C2 and C4. Notably, C1 exhibited the highest abundance of CD8+ T cells and the highest cytotoxicity score, indicating a strong anti-tumor immune response. C3 showed the highest abundance of cancer-associated fibroblasts (CAF) and endothelial cells, suggesting a tumor-promoting microenvironment. In addition, C1 and C3 had the highest abundance of myeloid dendritic cells and macrophages, which have complex roles in antigen presentation and tumor progression.

Further investigation into the expression of human leukocyte antigen (HLA) genes and immune checkpoint genes (ICGs) revealed significant upregulation of multiple HLA genes and ICGs in C1, including HLA-A, HLA-B, HLA-C, PDCD1, CD274, CTLA-4, HAVCR2, LAG3, CD80, CD86, TIGIT, etc. In contrast, although C3 had upregulated expression of HLA genes, its ICGs expression was weak ( Supplementary Figure S4B ).

These findings indicate that C1 is in an immunosuppressive state, potentially making it more responsive to immune checkpoint blockade (ICB) therapies by alleviating this inhibition and thereby inducing an immune response.

Considering that SNV can drive tumor progression by altering the structure of encoded proteins or regulating gene expression, therefore, we investigated the SNV landscape among clusters.

The SNV analysis showed that C1 exhibited high-frequency TP53 missense mutations (29.58%), low frequencies of TTN (15.49%) and CTNNB1 (15.49%) missense mutations; C2 harbored multiple high-frequency mutations, especially TTN (28.00% missense mutations and 13.00% silent mutations) and CTNNB1 (43.00% of missense mutations), with moderate TP53 mutations (13.00% of missense mutations); Notably, C3 had the lowest overall mutations among the four clusters; C4 was mainly characterized by high frequencies of TP53 (33.78% missense mutations) and TTN (25.68% missense mutations) mutations, while the frequency of CTNNB1 mutations was not significantly increased (27.03% missense mutations) ( Figure 5 ).

To evaluate the differences in drug sensitivity among clusters, we calculated the half-maximal inhibitory concentration values for 198 drugs. The results revealed distinct drug sensitivities for each cluster. C1 exhibited increased sensitivity to several drugs: gemcitabine, which interferes with DNA synthesis; cisplatin, which induces DNA cross-linking and damage; irinotecan, which inhibits topoisomerase I; and talazoparib, which inhibits poly ADP ribose polymerase (PARP). C2 was primarily sensitive to TAF1 inhibitor TAF1_5496_1732, which inhibits transcription; MCL-1 inhibitor AZD5991_1720, which promotes apoptosis; and mitochondrial inhibitor Dihydrorotenone_1827. C3 showed heightened sensitivity to nutlin-3a, which inhibits MDM2; selumetinib, which inhibits MEK1/2; mitoxantrone, which inhibits topoisomerase II, and gemcitabine. C4 was sensitive to Sepatronium bromide_1941, which inhibits survivin; Daporinad_1248, which inhibits nicotinamide phosphoribosyltransferase (NAMPT); and MK−1775_1179, which inhibits WEE1 and PLK1 in the cell cycle pathway ( Supplementary Figure S5A ; Supplementary Table S7 ).

We further investigated the correlation between drug sensitivity and TP53 mutations. The results showed that the expression of TP53 and MDM2 was higher in wild-type TP53 patients compared to those with TP53 mutations. Most patients in the C3 subtypes met the conditions of wild-type TP53 and high expression of MDM2, which could explain the sensitivity of nutlin-3a in C3 patients ( Supplementary Figures S5B, C ).

To gain deeper insights into clusters at the cell-subpopulation level, we performed single-cell analyses. Following quality control, dimensionality reduction, clustering, and cell annotation, we identified eight distinct cell subpopulations ( Supplementary Figure S6 ; Figures 6A, B ).

We evaluated the significantly upregulated marker gene sets of C1–C4 using the GSVA algorithm in single cells ( Figure 6C ). The results showed that the activity of the characteristic gene set of C1 was upregulated in various immune cell types, yet exhibited downregulated activity in malignant cells. In contrast, the activity of the characteristic gene set of C2 was significantly upregulated in malignant cells. The characteristic gene set of C3 showed pronounced activation in CAF, whereas the C4 characteristic gene set displayed activation across various cell types, predominantly in TEC, but showed diminished activity in B cells and CAF.

Cell-cell communication analysis revealed that tumor-associated macrophages (TAM) played a pivotal role in intercellular interactions within HCC ( Figure 6D ). The intercellular-interaction gene sets were used as reference gene sets for the GSEA enrichment analysis of C1–C4 based on bulk RNA-seq clustering ( Figure 6E ). The results showed upregulations of TAM-mediated and CAF-mediated intercellular communication in C1, a downregulation of CAF-mediated intercellular communication in C2, and an upregulation of intercellular communication mediated by CAF in C3 (P < 0.05). However, there was no significant enrichment of intercellular communication observed in C4.

To efficiently identify the molecular subtypes of patients in real-world settings, we developed a classifier using the random forest (RF) algorithm ( Figure 7A ).

We identified 521 cluster-specific marker genes using the “FindAllMarkers” function ( Supplementary Table S8 ). Based on these marker genes, the RF classifier achieved an overall prediction accuracy of 81.82% with a 95% confidence interval (CI) of 73.33%–88.53% in the TCGA internal validation cohort. The prediction accuracy for each subtype by the RF classifier is presented in Figure 7B .

In the ICGC external validation cohort, the molecular classifications were predicted using the RF classifier. Compared to the TCGA cohort, the predicted classifications in the ICGC cohort exhibited consistent prognostic characteristics (median OS: C3 > C2 > C1 > C4, P < 0.001, Figure 7C ). We also evaluated the characteristics of FSs ( Figure 7D ), immune cell infiltration ( Supplementary Figure S7A ), and immune genes ( Supplementary Figure S7B ) in the ICGC cohort. These validation results were consistent with the characteristics in the TCGA cohort, which further confirmed the robustness and reliability of our classifier.

Using our classifier, we further tested the hypothesis that C1 is more sensitive to anti-PD1 therapy in two external cohorts including 81 patients. Among them, the PRJEB34724 cohort enrolled 40 patients (6 responders and 34 non-responders), and the GSE202069 cohort enrolled 17 patients (8 responders and 9 non-responders).

The results showed that the immune cell infiltration and ICG expression profiles among the four subtypes were highly similar to the corresponding subtypes in the TCGA cohort ( Figure 8A ). Notably, we confirmed that patients in the C1 subgroup had a higher anti-PD1 response rate of 47.37% (9/19, P = 0.01), especially in the GSE202069 cohort with a response rate of 75% (6/8). Moreover, none of the C3-subtype responders were found in either of the two cohorts ( Figures 8B, C ).

To develop the prognostic nomogram, we identified 70 genes based on the intersection of cluster-specific marker genes and HCCDB prognostic genes. Then seven models were evaluated using the random survival forest (RSF) ( Supplementary Figures S8A–C ), the least absolute shrinkage and selection operator (LASSO) combining 10-fold cross-validation ( Supplementary Figures S8D, E ), and stepwise Cox regression. The details of covariables included in the seven models are shown in Supplementary Table S9 .

Using the 3-fold 1000-time-repeation cross-validation, the median concordance index (C-index) of the LASSO combined with stepwise Cox regression and tumor staging (LASSO-Stage) model achieved 0.770 in the training cohort and 0.810 in the TCGA inner validation cohort ( Supplementary Figure S8F ). In addition, we compared the time-dependent receiver operating characteristic curves and found that the LASSO-Stage model had higher values of area under curves at different time points ( Supplementary Figures S8G–J ). Therefore, the 9-variable LASSO-Stage model was relatively simple and accurate, and accordingly was used to develop the prognostic nomogram ( Figure 9A ). The C-index of the nomogram was 0.782 (95% CI, 0.709–0.855) in the external ICGC validation cohort. The calibration curves demonstrated good consistency between the nomogram predictions and actual outcomes ( Figures 9B–E ).

Decision curve analysis indicated that the 9-variable nomogram model was almost as effective as the 26-variable Step-Cox model. Patients could benefit from the nomogram when the risk thresholds range from 0.036 to 0.66 ( Supplementary Figures S9A, B ).

Using the nomogram, the risk scores were calculated for all patients, and the iterative stratification method was used to divide the training cohort into three stratifications (P < 0.001). Using the thresholds of the training cohort, patients in the validation cohorts could be divided into significant prognostic stratifications (P < 0.01) ( Supplementary Figures S10A–C ). The risk curves and survival status demonstrated that the risk scores of patients predicted by the nomogram were positively related to their death risk ( Supplementary Figures S10D–F ).

To facilitate classification and prognosis prediction, we developed an R package called “oncoClassSurv” (https://github.com/OliveryYL/oncoClassSurv), which has a corresponding desktop executable version (https://github.com/OliveryYL/oncoClassSurv_Expansion/) and a web-based Shiny calculator (https://oncomanager.shinyapps.io/oncoClassSurv/).

We identified four crucial prognostic genes through Venn analysis of multiple models ( Figure 10A ). Protein expression profiling and immune fluorescence indicated that FTCD was significantly downregulated in HCC (P < 0.001, Figures 10B, C ). Survival analysis revealed that patients with elevated FTCD expression had a significantly improved prognosis compared to those with diminished expression in the TCGA cohort (P < 0.001, Figure 10D ). FTCD expression was relatively upregulated in C2 and C3 with a favorable prognosis, while it was downregulated in C1 and C4 with a poorer prognosis ( Supplementary Figure S11 ).

We established FTCD-overexpressing Huh7 cell lines (Huh7-FTCD) and negative control (Huh7-NC) via lentiviral transduction ( Figure 10E ). Cell proliferation assays indicated that FTCD overexpression significantly inhibited the viability of Huh7 cells ( Figure 10F ). Colony formation assays demonstrated a reduced clonogenic ability in Huh7-FTCD cells ( Figure 10G ); Transwell assays revealed a significantly decreased migration ability in Huh7-FTCD cells ( Figure 10H ). Apoptosis assays showed no significant difference between Huh7-FTCD and Huh7-NC cells ( Figure 10I ), while the Huh7-FTCD cells exhibited G2/M cell cycle arrest compared to Huh7-NC cells ( Figure 10J ). Subcutaneous tumorigenicity assay in nude mice indicated that tumor progression originating from the Huh7-FTCD cells was significantly slower than that from the Huh7-NC cells ( Figure 10K ).