A total of 15 independent public datasets were collected and processed in this study, including single-cell RNA-seq cohort, high-throughput RNA-seq cohorts, microarray cohorts, and clinical treatment cohorts. The single-cell RNA-seq cohort GSE140228 was downloaded from Gene Expression Omnibus (GEO), which deciphered the immune landscape and dynamics of HCC. The TCGA-LIHC cohort (n = 369), International Cancer Genome Consortium (ICGC)-LIRI cohort (n = 232), and GSE14520 cohort (n = 221) included gene expression, and corresponding complete clinical information were retrieved from TCGA, ICGC, and GEO, respectively. In the TCGA-LIHC cohort, somatic mutation data and copy number variation (CNV) information were accessed from the online portal cBioPortal. In addition, the E-TABM-36 cohort (n = 60) was generated from ArrayExpress database. Other microarray cohorts were also collected from the GEO database, encompassing GSE25097 (n = 268), GSE76427 (n = 115), GSE116174 (n = 64), GSE144269 (n = 68), GSE104580 (n = 147), and GSE109211 (n = 67). Among these, GSE104580 and GSE109211 contained TACE treatment and sorafenib therapy information, respectively. Additionally, we enrolled four eligible immunotherapy cohorts with 98 non-responders and 41 responders, including GSE35640 (n = 56), GSE91061 (n = 39), GSE100797 (n = 21), and Nathanon (n = 23) cohorts. The details of all retrieved cohorts are shown in Supplementary Table S1.

The RNA-seq raw count data were converted to transcripts per million format and further log-2 transformed. The expression profiles from GEO and ArrayExpress databases were processed and normalized by affy and lumi packages based on different platforms. According to the RECIST v1.1 standard, patients with complete response/partial response and patients with stable disease/progressive disease were deemed as responders and non-responders, respectively. Patients who were not evaluable were excluded.

The Seurat (v4.0.6) package was utilized to process data for further dimension reduction and cell clustering analysis (19). Single-cell gene expression profiles were filtered to exclude cells that had either over 10% mitochondria genes or fewer than 200 transcripts/cell. PCA linear dimensional reduction and clustering visualization were performed using RunPCA function and RunTSNE function implemented in Seurat. The SingleR package was applied to annotate distinct cell clustering, and then unique marker genes were identified via the FindAllMarkers function of Seurat. The irGSEA package with UCell method was used to accomplish single-cell gene set enrichment analysis. The pseudo-time trajectory analysis of single cells was conducted by Monocle 2 package (20). Cellular trajectory ordered in pseudo-time was presented with multiple branches, and genes along the trajectory were enrolled in a subsequent analysis.

Nonnegative matrix factorization (NMF) algorithm executed in the NMF package was performed to identify molecular clustering by factorizing matrix and running iterations (21). Using univariate Cox regression, the trajectory genes were screened, and prognosis-associated candidate genes were generated for a better clinical application. Based on the nonnegative matrix of these genes, consensus clustering was deciphered with the following criteria in the NMF package: possible factorization ranks = 2–9, number of iterations = 100, and method = “lee”. Cophenetic coefficient was employed to determine optimal rank, and silhouette statistic was used to quantify the robustness of clustering patterns. Usually, when the value of cophenetic correlation coefficient starts decreasing, it is deemed as optimal factorization rank (21). The magnitude of silhouette coefficient was linked with the similarity of a sample to its own cluster, a higher silhouette value, and a better one matched to its own cluster (22).

The weighted gene co-expression network analysis (WGCNA) aims to explore and reveal the correlations between gene modules and phenotypes (23). The characteristic genes of distinct clusters were identified using the WGCNA package. After excluding the outlier samples, gene co-expression network was constructed based on the top 5,000 genes and further transformed into a scale-free network via selecting an appropriate soft threshold β. Subsequently, the topological overlap matrix (TOM) describing the overlap of network neighbors and 1-TOM representing gene dissimilarity were generated by the weighted adjacency matrix. Eventually, gene modules with various colors were identified by dynamic tree algorithm. Based on the relationship between the module eigenvalue and phenotypes, three modules with highest correlation were filtered, and the characteristic genes were acquired for subsequent analysis.

The nearest template prediction (NTP) is a flexible approach that evaluates class prediction confidence for a single patient (24). To further assess the reliability and stability of clusters, the NTP algorithm implemented in the CMScaller package was utilized to validate by multiple cohorts from inconsistent platforms. The signature gene list used in NTP was derived from modules’ characteristic genes and differentially expressed genes.

To explore the specific biological characteristics of distinct clusters, gene set enrichment analysis (GSEA) algorithm was performed, which displayed pathway activities by gene rank information. Differentially expressed genes (DEGs) were identified via the limma package and then ordered according to descending log₂ fold change value. The clusterProfiler package was used to exhibit the GSEA analysis, and Benjamin–Hochberg-corrected adjusted P-value <0.05 was regarded as statistically significant. The gene set variation analysis (GSVA) was broadly utilized in pathway activity assessment (13). Based on 50 Hallmark gene sets, the GSVA package was applied to further evaluate and elucidate different biological characteristics among clusters.

The landscape of genomic variations was depicted by the mutation and CNV data. The maftools package was used to display somatic variants among distinct clusters, including single-nucleotide polymorphism (SNP), insertion and deletion (INDEL), tumor mutation burden (TMB), and mutation frequency (25). Generally, frequently mutated genes (FMGs) that had top 20 mutation frequency were considered the main driver genes for malignant tumor (26). We also dissected the CNV among clusters and further exhibited frequently AMP or HOMDEL genes, which possessed the top 10 genes with amplification or deletion.

To decode the landscape of immune cell infiltration, immune gene sets were obtained from a previous study which stored 28 immune cell subgroups (27) (Supplementary Table S2). The single-sample gene set enrichment analysis (ssGSEA) algorithm was used to compute the relative infiltration abundance of 28 immune cells. The expression of 27 immune checkpoints was evaluated to further depict the tumor immune microenvironment, encompassing the B7-CD28 superfamily, TNF superfamily, and other molecules (15) (Supplementary Table S3). The capacity of antigen presentation was estimated based on nine human leukocyte antigen (HLA) molecule expression profiles (16). Two prevalent approaches were employed to assess the immunotherapeutic efficacy among distinct clusters, including T cell inflammatory signature (TIS) and unsupervised subclass mapping (Submap). TIS, which contained 18 inflammatory genes, was scored by ssGSEA algorithm, which gave a higher score to indicate a better response to PD-1 blockade (28). The Submap was utilized to measure the expression profile similarity between HCC patients and immunotherapeutic patients, which was consistent with the similarity of clinical responses (29). Four immunotherapy cohorts were applied to further reveal the immunotherapy significances of different clusters. Receiver operating characteristic curve (ROC) was executed to estimate the accuracy of the immunotherapeutic prediction.

The pRRophetic package encompassing linear ridge regression model was applicable to predict drug response based on gene expression data (30). Using pRRophetic package, we calculated the half-maximal inhibitory concentration (IC50) of clinical tissues. Potential therapeutic agents were screened when the IC50 value was lowest among distinct clusters. The Connectivity Map (CMap) was a systematic approach for searching potential therapeutic compounds based on the similarity of gene expression profile (31). We executed CMap database to identify potential therapeutic compounds and target pathways. DEGs were firstly screened by the limma package, and the expression similarity was compared with database signatures, and then enrichment score-quantified therapeutic value was generated. Additionally, the TACE and sorafenib treatment-associated cohorts were also enrolled to assess the clinical efficacy among distinct clusters.

The Kaplan–Meier and Cox regression analyses were conducted by the survival package. The log-rank test was applied to compare the survival statistics of categorical variables. Multivariate Cox regression analysis was utilized to calculate the hazard ratio and verify the independent significance of multiple traits. Using Pearson’s correlation analysis, the correlation between two continuous variables was evaluated. Kruskal–Wallis test was performed to compare the difference among three groups. The pROC package was carried out to draw the ROC for predicting binary categorical variables. All data cleaning, statistical analysis, and visualization were conducted in R v4.1 software. All statistical tests were two-sided. P <0.05 was considered statistically significant.

To depict the landscape and dynamics of CD45+ immune cells, single-cell transcriptomes were performed with subcluster analysis and visualized by t-distributed stochastic neighbor embedding (t-SNE) approach. A total of 3,817 immune cells originated from tumor and normal tissue were classified into 17 clusters (Figure 1A). Based on the gene signature of distinct subclusters, these cells were mainly composed of five immune cell clusters, including B cells, dendritic cells, monocyte cells, natural killer cells, and T cells (Figure 1B). To describe the source of immune cells, 2,212 tumor-derived cells and 1,605 non-tumor-derived cells were clustered separately (Figure 1C). Single-cell enrichment analyses exhibited that all cells were involved in immune inflammation pathway, especially dendritic cells, monocyte cells, and T cells (Figure 1D). We further measured the percentage composition and extracted the top 5 markers of each immune cells cluster (Figures 1E, F). In addition, the accuracy of B cell populations was verified by analyzing the expression of specific markers: CD79A, IGHG3, and IGLC2 (Supplementary Figures S1A–F). A concern was that T cells from tumor tissue displayed the most prevalent cell cluster. It was well known that T cells had anti-tumor immunity ability and power in directly killing cells in tumor progression. Thus, T cells were further explored and decoded at single-cell level.

The tumor-derived T cells were employed to reveal the evolution of T cells using dimensional reduction, unsupervised clustering, and trajectory analysis. All T cells were clustered again and divided into six cell subpopulations (Figure 2A and Supplementary Figure S2A). A previous study has elucidated that CD3D and CD3E were shared gene markers of CD4+ and CD8+ T cells (32). Thus, our study suggested that the T cells retrieved in this study mainly consisted of CD4+ and CD8+ T cells (Figures 2B, C and Supplementary Figures S2B, C). The specific marker gene CD4 was expressed in clusters 1, 2, 3, and 5, which indicated that these clusters represented CD4+ T cells (Figure 2D and Supplementary Figure S2D). We also noticed that cluster 0 and cluster 4 were enriched for CD8+ T cell markers, such as CD8A and CD8B, confirming the identity of CD8+ T cells (Figures 2E, F and Supplementary Figures S2E, F). These specific cell subpopulations were visualized with two-dimensional distributions by the t-SNE method (Figure 2A). As the tumor progressed, the cell status in TME was dynamic rather than immobile, and the dynamic process was depicted by the Monocle algorithm. For CD4+ T cells, there was a gradual evolution process from CD4+ T clusters 3 and 4 to CD4+ T clusters 1 and 2 (Supplementary Figures S2G, H). The immune checkpoints were further investigated in CD8+ T cells (32). Strikingly, most inhibitory checkpoints represented T cell exhaustion, such as HAVCR2 (TIM3), LAG3, TIGIT, PDCD1 (PD-1), and CTLA-4, which were significantly upregulated in CD8+ T cluster 2 (Figure 2G). The pseudo-time and trajectory analysis indicated that CD8+ T cells tended to be exhausted with tumor progression, which might be linked with poor prognosis (Figures 2H, I). A previous study had also demonstrated that T cell exclusion is common in TME and associated with immune privilege (33). Thus, these genes along the trajectory might play an important role in TME and in the clinical outcomes of HCC patients (Supplementary Table S4).

The trajectory genes were employed to extract prognosis-associated candidate genes. Based on the expression of these genes, the NMF approach was utilized to decipher heterogeneous molecular clusters. As illustrated in Figure 3A, the optimal cluster option was three due to the cophenetic coefficient that started to rapidly decline. The consensus matrices also suggested that three clusters had optimal stratification (Figure 3B). The silhouette statistic was used to assess the stability of molecular clusters, and the samples were further detected by silhouette width (34). Therefore, samples with a positive silhouette width were divided into three stable and robust clusters (Figure 3C). To facilitate the clinical application, the prognostic significance of clusters was further explored. C2 exhibited poor overall survival (OS) and recurrence-free survival, whereas C3 presented a favorable prognosis (P <0.05) (Figures 3D, E).

The characteristic genes of three clusters were identified using WGCNA package. First of all, the outlier samples were removed, and then the remaining samples were clustered (Supplementary Figure S3A). When β was set to 6, the no-scale R ² was 0.9, developing a scale-free network (Figure 4A). As shown in Supplementary Figure S3B, the TOM network was displayed via a heat map. Subsequently, a total of 12 co-expression modules were obtained by dynamic tree cutting, and the eigengene adjacency of various modules was depicted via a heat map (Supplementary Figures S3C, D). Furthermore, the module–trait relationships were exhibited to measure the correlations between the modules and the three clusters. The turquoise, blue, and purple module presented the strongest correlation with C1, C2, and C3, respectively (Figure 4B). The correlation values between gene significance and module membership indicated that the construction of the gene modules was robust (Figures 4C–E). The genes in each module were defined as characteristic genes and are shown in Supplementary Table S5. To decode the specific biological characteristics of three clusters, we performed GSEA analysis using gene sets from Gene Ontology and Kyoto Encyclopedia of Genes and Genomes. C1 was mainly associated with metabolism pathways, including fatty acid oxidation and bile secretion. There are tight links between tumor immune microenvironment and metabolism activity such as fatty metabolism and bile acid metabolism in liver. The important source of energy is generated from the elevated lipid metabolism activity, which is the power and key regulator for immune cells and tumor cells. Meanwhile, elevated lipid metabolism activity impacts the inflammatory pathway in the tumor microenvironment, even resulting in immune escape (35). C2 possessed conspicuous enrichment in proliferation pathways such as cell cycle and nuclear division. C3 was obviously enriched in immune pathways encompassing the adaptive immune response and chemokine signaling pathway (Figures 4F, G). Therefore, we characterized C1 as lipid metabolic HCC and C2 as cell proliferative HCC, whereas C3 was defined as immune inflammatory HCC.

Based on signature gene expression, NTP analysis was performed to assess the reliability and stability of the three clusters. The signature genes were generated from the overlaps between characteristic genes and upregulated DEGs (24). Using the NTP method, a total of eight cohorts obtained from distinct platforms were executed to measure and evaluate the prediction confidence for each patient, including GSE14520 (Figure 5A), ICGC-LIRI (Figure 5C), E-TABM-36, GSE76427, GSE25097, GSE104580, GSE116174, and GSE144269 (Supplementary Figures S4A–F). Consistent with a previous study, patients with false discovery rate <0.05 were detected for a subsequent analysis (36). In GSE14520 and ICGC-LIRI cohorts, Kaplan–Meier and multivariate Cox regression analyses were utilized to further elucidate the prognostic implications of the three clusters. The Kaplan–Meier analysis suggested that C2 still possessed the most unfavorable OS, while C3 presented the most favorable OS (P <0.05), which was coincident with previous results (Figures 5B, D). Multivariate Cox regression indicated that C2 was an independent prognostic indicator in the TCGA-LIHC, GSE14520, and ICGC-LIRI cohorts (Figures 5E–G). In addition, the proportions of the three clusters were displayed among distinct cohorts, which showed a high similarity (Figure 5H). Overall, the three clusters had heterogeneous clinical outcomes and were reproducible and robust in HCC.

We further depicted the landscape of genomic variations among the three heterogeneous clusters. As illustrated in Figure 6A, the mutation frequency of the top 20 FMGs was exhibited, and the overview of SNP, INDEL, and TMB was also displayed. To increase the understanding of somatic mutation, we compared the mutational differences of 20 FMGs among three clusters (Figure 6B). Strikingly, there were three universal FMGs exhibited in all molecular clusters, encompassing TP53, CTNNB1, and TTN, indicating that these FMGs might play a key role in tumorigenesis or tumor progression (Figure 6B). Among the three FMGs, the CTNNB1 and TTN mutations were pronounced in C1, while the TP53 mutation was dominant in C2. A previous study had elucidated that the gradual accumulation of gene mutations was prone to tumorigenesis (37). C3 possessed the lowest gene mutation relative to others, implying better clinical outcomes. We further explored the summary of CNV and depicted the top 10 frequent AMP and HOMDEL genes among the three clusters (Figure 6C). Interestingly, C2 was characterized by a higher CNV compared to the other clusters, and CSMD1 gene had the most conspicuous CNV loss (Figure 6C). In line with previous research, the loss of a putative tumor suppressor gene, CSMD1, might be the driven event in HCC progression (38, 39). Taken together, patients in C2 conveyed prominent genomic variations, indicating a highly genomic instability.

The 50 Hallmark gene sets were broadly performed in cancer-related research (14). We further revealed the potential carcinogenic characteristics of the three clusters using GSVA algorithm. Consistent with the above-mentioned results, C1 was characterized by lipid metabolic pathways, and C2 was mainly associated with cell proliferative pathways, while C3 was enriched in immune inflammatory pathways (Supplementary Figure S5A). These results indicated that patients in C3 might obtain better immunotherapeutic efficacy. Therefore, the landscape of immune cell infiltration and immune checkpoint expression was delineated to decipher the underlying mechanism. Compared to other clusters, C3 exhibited more infiltration abundance of immune cells (Figure 7A). Patients in C3 tended to be the “immune-hot” subtype, which stored massive immune cells in TME, including CD4+ T cell, CD8+ T cell, activated dendritic cell, nature killer cell, and so on (P < 0.05) (Supplementary Figure S6A). Among the three clusters, C3 also displayed the highest expression of immune checkpoints, such as CD274 (PD-L1), CTLA-4, and LAG3, which suggested that C3 might be more sensitive to immune checkpoint inhibitors (ICI) therapy (Figure 7B). Furthermore, the expression of HLA molecules was conspicuously higher in C3, which proved that patients in C3 possessed a strong power of antigen presentation (Figure 7C). Overall, precision immunotherapy might be applicable to patients in C3.

To yield deep insights on immunotherapy, the TIS and Submap methods were applied to assess the clinical efficacy among distinct clusters. As expected, C3 had the highest TIS score, hinting immune activation and elevated response to ICIs (P < 0.0001) (Figure 7D). Based on the Submap algorithm, patients in C3 had a superior similarity of expression patterns with patients responding to PD-L1 inhibitor (Bonferroni corrected P <0.01), indicating more benefits from anti-PD-L1 treatment (Figure 7E). Subsequently, the characteristic genes of C3 were extracted to serve as the measured score, named as cluster-associated immunotherapy score (CAIS). The patients were classified into high and low groups based on the constant ratio (3 vs. 7). Four immunotherapy cohorts contained 98 non-responders, and 41 responders were retrieved to evaluate the clinical applications of CAIS. Notably, patients in the high group exhibited a superior response to immunotherapy in all cohorts (Figure 7F). The area under the curve values of CAIS for predicting the accuracy of immunotherapeutic efficacy were 0.721, 0.709, 0.719, and 0.737 in GSE100797, GSE35640, GSE91061, and Nathanon cohorts, respectively (Figure 7G). Collectively, the above-mentioned data suggested that CAIS was a robust and promising immunotherapy indicator. Patients in C3 were recommended to be taken into consideration in using immunotherapy.

Two other clinical treatment cohorts, GSE140580 and GSE109211, were collected to seek a high-potency treatment strategy and facilitate the clinical benefits for HCC patients. The TACE and sorafenib therapies were widely applied to clinical practice. As displayed in Supplementary Figure S5B, both C1 and C3 were characterized with early AJCC stage, which is predominantly associated with a superior prognosis. According to the updated clinical guideline, the TACE treatment was recommended for patients with primary stage (7). In line with that, our results also demonstrated that patients in C1 and C3 had a desirable efficacy for TACE treatment (Figure 8A). The therapy sensitivity of sorafenib was further estimated among the three clusters, and it was substantiated that patients in C2 could obtain more clinical benefits (Figure 8B). To identify latent therapeutic drugs, the pRRophetic package was used to evaluate the sensitivity of numerous agents, which was quantified by half-maximal inhibitory concentration (IC50). In the TCGA-LIHC cohort, patients in C2 also displayed superior response to nilotinib and bosutinib relative to other clusters (Figures 8C, D). Furthermore, patients in C1 might be more sensitive to axitinib for the lower IC50 value (Figure 8E). Obatoclax, docetaxel, and cisplatin were potential therapeutic agents for patients in C3 (Supplementary Figures S6B–D). These drugs might be developed into promising therapeutic agents for distinctly classified HCC patients. The CMap database contained a massive gene expression profile, which could assess the relationships among gene expression, phenotype, and drugs, and was further utilized to identify other potential compounds and shed light on the mode of action. We exhibited 20 compounds that possessed personalized therapeutic potential for three clusters (Figure 8F). The targeted pathways of these compounds were depicted, which could be used to develop numerous drugs (Figure 8G). According to patients in distinct clusters, applying individualized therapy patterns would produce curative clinical efficacy.