Integrating single-cell and bulk RNA sequencing data to characterize the heterogeneity of glycan-lipid metabolism polarization in hepatocellular carcinoma

Abstract

Background

Hepatocellular carcinoma (HCC) is high heterogeneity and remains an unmet medical challenge, but their metabolic heterogeneity has not been fully uncovered and required clinical applicable translational strategies.

Methods

By analyzing the RNA sequencing data in the in-house cohort and public HCC cohorts, we identified a metabolic subtype of HCC associated with multi-omics features and prognosis. Multi-omics alterations and clinicopathological information between different subtypes were analyzed. Gene signature, radiomics, contrast-enhanced ultrasound (CEUS), serum biomarkers were tested as potential surrogate methods for high throughput technology-based subtyping. Single-cell RNA sequencing analyses were employed to evaluate the immune characteristics changes between subtypes.

Results

By utilizing metabolic-related pathways, we identified two heterogeneous metabolic HCC subtypes, glycan-HCC and lipid-HCC, with distinct multi-omics features and prognosis. Kaplan–Meier and restricted mean survival time analyses revealed worse overall survival in glycan-HCCs. And glycan-HCCs were characterized with high genomic instability, proliferation-related pathways activation and exhausted immune microenvironment. Furthermore, we developed gene signatures, radiomics, CEUS and serum biomarkers for subtypes determination, which showed substantial agreement with high-throughput-based classification. Single-cell RNA-seq showed glycan-HCCs were associated with multifaceted immune distortion, including exhaustion of T cells and enriched SPP1 + macrophages.

Conclusion

Collectively, our analysis demonstrated the metabolic heterogeneity of HCCs and enabled the development of clinical translation strategies, thus promoting understanding and clinical applications about HCC metabolism heterogeneity.

Introduction

With 865,269 cases diagnosed and 757,948 deaths in 2022, liver cancer is ranked as the sixth most common tumor and the third leading cause of cancer-related death [1]. Hepatocellular carcinoma (HCC) accounts for approximately 90% of liver cancer cases [2]. HCC exhibits a high degree of complexity and heterogeneity in molecular phenotypes. Recently, multi-omics data such as genomics, transcriptomics, proteomics and epigenetics data, have been used to determine molecular subtypes of HCC and lead to deep insights into HCC [3,4,5]. These studies have highlighted the value of molecular features for HCC patients’ risk stratification and individualized treatment decision.

Metabolic reprogramming is a hallmark of cancer and its vulnerabilities has provided a challenge to developing therapies [6]. The liver is the largest organ in the body responsible for metabolic control, and metabolic disorders are closely associated with the occurrence of liver cancer [7]. The incidence of HCC caused by metabolic dysfunction-related factors is increasing [8]. Notably, glycan and lipid metabolism play pivotal roles in HCC cancer cell signaling, immune evasion, and metastasis [9, 10]. In contrast to other molecular classification systems, the metabolism-based classification system can provide more targeted guidance for metabolic therapy in HCC. Recent studies utilizing multi-omics data have elucidated the metabolic heterogeneity of HCCs and established metabolism-related subtypes of HCC. For example, Gao et al. [3] provided proteogenomic landscape of HCC and found three proteomic subgroups have distinct metabolic status. Li et al. [11] determined multi-omics classification of HCC based on the fatty acid degradation pathway and provided personalized therapy stranguries. However, most study for HCC molecular stratification is mainly based on high-throughput technology, it is limited by high costs, complex technical process, and potential batch effects. It is imperative to uncover its clinical value and boost clinical application.

Given the critical roles of glycan and lipid metabolism in HCC, this study aims to elucidate their roles in HCC tumor heterogeneity and provide multi-omics alterations insights into glycan-lipid metabolism polarization. Furthermore, our results proposed explored gene signature, medical images and serum biomarkers for subtypes alternative techniques. Single-cell RNA-seq data further provided immune landscape influenced and altered by metabolic heterogeneity. This exploration into the metabolic heterogeneity of HCC not only refines our understanding of HCC intrinsic properties but also provides clinically practical methods for subtypes determination.

Methods

In-house HCC cohort acquisition

The study was approved by the Research Ethics Committee of the First Affiliated Hospital of Guangxi Medical University, and written informed consent was obtained from each patient. Tumor tissues were surgically, or biopsy collected from 110 patients pathologically diagnosed with primary HCC at the First Affiliated Hospital of Guangxi Medical University and immediately stored in liquid nitrogen. Detailed clinical information is presented in Table 1. Total RNA was extracted from each sample using TRIzol reagent. High-quality RNA was used to construct sequencing libraries with the Illumina platform. The library preparation process involved poly-A selection or ribosomal RNA depletion for mRNA enrichment, followed by reverse transcription to synthesize cDNA, end repair, adaptor ligation, and PCR amplification. Sequencing was carried out on the NovaSeq 6000 platform, performing 150 bp paired-end reads with a target depth of approximately 30 million reads per sample. The raw sequencing data were initially quality checked via FastQC software. Adapters and low-quality reads were trimmed using Trimmomatic. Clean reads were aligned to the human reference genome (GRCh38) using HISAT2, with alignment rates and coverage assessed. Gene expression quantification was performed using featureCounts, and the read numbers per gene were normalized to the FPKM values.

Table 1 Clinicopathological characteristics of included in in-house cohort

Publicly available data collection and preprocess

We also collected three publicly available RNA-seq from The Cancer Genome Atlas Program (TCGA) (project ID: TCGA-LIHC) [4], International Cancer Genome Consortium (ICGC) (project ID: LIRI-JP) and Clinical Proteomic Tumor Analysis Consortium (CPTAC) (project ID: CHCC-HBV) [3] database to validate our subtyping strategies. For TCGA and ICGC datasets, we removed recurrent HCC tissues in this study.

Metabolic pathways score and consensus clustering

Metabolic pathways were downloaded from scMetabolism software, which collected 85 metabolism related pathways from KEGG (Kyoto Encyclopedia of Genes and Genomes) database (Supplementary Table S1) [12]. There pathways were clustered into 11 major categories based upon KEGG classifications, included Glycan biosynthesis and metabolism, Nucleotide metabolism, Amino acid metabolism, Carbohydrate metabolism, Lipid metabolism, Metabolism of cofactors and vitamins, Metabolism of other amino acids, Energy metabolism, Biosynthesis of other secondary metabolites, Xenobiotics biodegradation and metabolism, and Metabolism of terpenoids and polyketides. Single-sample gene set enrichment analysis (ssGSEA) method in gene set variation analysis (GSVA) package [13] was utilized to calculate the enrichment score of each metabolic pathway based on transcriptomic data.

Considering the sample size of the four cohorts, we utilized TCGA-LIHC cohort as the discovery cohort. TCGA cohort included 374 tumor samples and 50 non-tumor samples. Then, 85 metabolic pathways’ enrichment score of 371 primary HCC were submitted to consensus cluster. We employed unsupervised clustering analysis from R package “ConsensusClusterPlus” with parameters “clusterAlg = “pam”, reps = 1000, pItem = 0.8, pFeature = 1 and distance = euclidean” [14]. The optimal number of clustering was determined by the average pairwise consensus matrix and the proportion of ambiguous clustering (PAC). Then, for three validation cohort, included ICGC (n = 231), CPTAC (n = 159), and in-house cohorts (n = 110), we selected top 10 most up-regulated significant pathways for each cluster to achieve subtypes determination. Parameters consistent with those in TCGA cohort are used for clustering and clusters were determined based on pathways’ activation status.

From Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/index.jsp), KEGG pathways annotation was downloaded. These annotations were used for functional enrichment analysis. To perform GSEA, the differentially expressed genes were submitted to the clusterProfiler software to identification of statistically enriched pathways by comparing the distribution of DEGs against the KEGG pathway gene sets [15]. Four genomic instability related features, included aneuploidy score, copy-number alterations (CNA), homologous recombination deficiency (HRD) score and mutation burden were obtained from previous publication [16]. Microenvironment Cell Populations-counter (MCP-counter) method was used to calculated immune infiltrates in tumor samples [17]. We used average expression of 5 exhausted markers (CTLA4, HAVCR2, LAG3, PDCD1, and TIGIT) to define the exhausted score for T cells.

Contrast-enhanced ultrasound (CEUS) Liver Imaging Reporting and Data System (LI-RADS)

After the target lesion was identified by two-dimensional ultrasound, 1.5 mL of SonoVue contrast agent (Bracco Imaging Milan, Italy) was injected into the cubitus vein, followed by 5.0 mL of saline. Dynamic images of CEUS observed for either 3–5 min or until the microbubbles disappeared. Based on the 2017 edition of CEUS LI-RADS, lesions were divided into different categories of LR-1, LR-2, LR-3, LR-4, LR-5 and LR-M. Case inclusion criteria: (1) HCC included in our RNA-seq analysis; (2) Received CEUS examination within 2 weeks before surgery; (3) HCC high-risk patients, namely cirrhosis, chronic hepatitis B virus infection, current or past HCC. Exclusion criteria: (1) The sonography images are of poor quality; (2) Patients who had received preoperative transcatheter arterial chemoembolization, chemotherapy, radiotherapy, etc. The full version of CEUS LI-RADS can be found on the ACR website (https://www.acr.org/Clinical-Resources/Reporting-and-Data-Systems/LI-RADS). Radiologists were unaware of the patients’ clinical information and pathological results in evaluating.

Radiomics analysis

Available Computed Tomography (CT) images in in-house (n = 110) and TCGA cohorts (n = 40) were include for radiomics analysis. Three-dimensional region of interests (ROIs) of tumors were manually delineated by radiologists. Then, 107 radiomics features were extracted from ROIs by using Pyradiomics software, included 18 first-order, 14 shape and 75 textural features [18]. Wilcoxon tests were used to determine radiomics features that were differentially expressed (p < 0.05) between glycan-HCCs and lipid-HCCs. Then, these features were submitted to least absolute shrinkage and selection operator (LASSO) algorithm for radiomics model development. By imposing a penalty on the absolute values of the regression coefficients, LASSO effectively shrinks the coefficients of less relevant features to zero, thereby retaining only the most informative features and mitigating the risk of overfitting. The tuning parameter (λ) in the LASSO model was selected using ten-fold cross-validation, with the area under the curve (AUC) as the performance metric. The optimal λ, denoted as lambda.min to ensure the model achieved the highest discriminative ability. Finally, a radiomics score was calculated for each patient via a linear combination of selected features that were weighted by their respective coefficients.

scRNA-seq data analysis

Single-cell RANs-seq data and corresponding cell type annotation which deposited in the Gene Expression Omnibus (accession ID: GSE151530), including 46 HCC and intrahepatic cholangiocarcinoma biopsies from 37 patients, were downloaded for analysis. We included pretreatment samples from HCC patients to explore tumor immune microenvironment (TIME) differences between glycan-HCC and lipid-HCC. And samples with a tumor cell proportion greater than zero were finally included in the analysis (n = 14). Details of single-cell RNA sequencing (scRNA-seq) processes have been reported [19, 20]: Briefly, core needle biopsies or resected tumor tissues were dissociated using the Miltenyi Tumor Dissociation Kit and processed into single-cell suspensions. Libraries were prepared with the 10× Genomics, and analyzed using Cell Ranger. Malignant and non-malignant cells were distinguished by inferring copy number variations (CNVs). And cell annotations were based on metadata file provided by dataset.

Then, ssGSEA was used to calculate metabolic pathway enrichment score of each cell. For each sample, the average enrichment score of all malignant tumor cells was calculated and used as the representative metabolic pathway enrichment score for that sample. These scores were then subjected to supervised clustering to classify samples into distinct metabolic subtypes, specifically glycan-HCC and lipid-HCC [21]. By evaluating the ligands, receptors, and their respective co-factors to determine communication information between glycan and lipid HCCs.

Statistics analysis

All statistics analysis was performed by using R software (version 4.3.0). Genes and pathways differentially expressed between two clusters were used Wilcoxon test. Gene mutation status between different groups were compared by using the chi-square test. For survival analysis, Survival curves were generated using the Kaplan–Meier (K–M) method to estimate the OS probabilities. Differences between groups were compared using the log-rank test, which assesses the significance of differences in survival distributions. Hazard ratios (HRs) and 95% confidence intervals (CIs) were reported to quantify the association between subgroups and outcomes. Restricted mean survival time (RMST) analysis was also performed using the “survRM2”. The cut-off values were established based on clinical guidelines and statistical methods. Specifically, for AFP, the cut-off value was determined using established clinical thresholds (≥ 400 ng/ml). The optimal cut-off value for radiomics was determined using Youden’s Index. A two-sided p value < 0.05 was considered the statistical significance threshold and the Bonferroni–Holm correction to adjust for multiple testing.

Results

Overview of study design

In this study, we sought to characterized and boost clinical translational of HCC metabolic-pathways-based subtypes. To achieve this goal, this study mainly included four phases (Fig. 1). In the first phase, we utilized in-house and three public RNA-seq datasets to develop and validate two metabolic-pathways based subtypes. In the second phase, we determined the prognostic values, and multi-omics features between two metabolic subtypes. In the third phase, gene signature algorithm, CEUS LI-RADS, radiomics technology and serum AFP were determined as surrogated methods for high-throughput technology developed subtyping algorithm. In the last phase, we utilized single-cell RNA-seq analysis to determine relationships between metabolic subtypes and TIME.

Fig. 1

Schematic overview of the study design

Transcriptome profiles defined metabolism polarization

Considering the sample size of the four cohorts included in our study, TCGA was used as discovery cohort. In TCGA cohort, unsupervised cluster showed distinct metabolism pattern in HCC (n = 374) and non-tumor (n = 50) tissues (Supplementary Figure S1A). UMAP also demonstrated that the enrichment score of metabolic-pathways were significantly altered between HCC and non-tumor tissues (Supplementary Figure S1B). Wilcoxon test revealed that “Glycan biosynthesis and metabolism” and “Nucleotide metabolism” were the two major pathways up-regulated in HCC tumor tissues. And other metabolic-pathways, included amino acid, lipid and carbohydrate metabolism were significantly down-regulated in tumor tissues when compared with non-tumor tissues (Supplementary Table S2). Then, enrichment score of the 85 metabolic-pathways in 371 primary HCC samples were submitted to consensus clustering for subtypes identification. We found the highest average consensus score at k = 2 (Supplementary Figure S1C). The PAC score reached its minimum at k = 2, demonstrating maximal clustering stability, and exhibited a significant increase at k = 3, thereby supporting k = 2 as the optimal clustering solution (Supplementary Figure S1D). And clear block-like structures in the consensus matrix, indicating robust cluster separation at K = 2 (Supplementary Figure S1E–G). Therefore, we finally determined the cluster number as 2.

Two subtypes showed distinct metabolic-pathways alteration status (Supplementary Figure S2; Supplementary Table S3). Cluster1 was designated as the glycan-HCCs subtype, considering it was characterized by remarkable upregulation of glycan metabolism pathways, including “Mannose type O-glycan biosynthesis”, “Glycosaminoglycan biosynthesis—heparan sulfate/heparin” and so on. Cluster2 was characterized as the relative upregulation of most metabolism-related pathways. Notably, lipid metabolism pathways, such as “Primary bile acid biosynthesis”, “Steroid hormone biosynthesis” and “Fatty acid degradation”, were the main components of the first 10 up-regulated pathways. And cluster2 was determined as lipid-HCCs.

Supervised analyses were performed in ICGC, CPTAC and in-house based on top 10 enriched pathways for glycan- and lipid-subtypes respectively (Fig. 2a). ICGC (Fig. 2b), CPTAC (Fig. 2c) and in-house (Fig. 2d) cohorts were separated into two subgroups. Then, according to the enrichment score of the 20 metabolic pathways in the two subtypes, they were determined as glycan or lipid subtypes in ICGC (Fig. 2e), CPTAC (Fig. 2f) and in-house (Fig. 2g) cohort respectively.

Fig. 2

Determination of metabolic subtypes in three external cohorts. a Top 10 most significantly upregulated pathways in glycan-HCCs and lipid-HCCs subtypes. b–d Consensus matrix heatmaps from supervised clustering analyses of ICGC (b), CPTAC (c), and in-house cohorts (d), generated using the ConsensusClusterPlus R package. e–g Subtype stratification based on pathway enrichment scores in ICGC (e), CPTAC (f), and in-house cohorts (g). Heatmap color gradients represent normalized enrichment scores of subtype-specific pathways. Heatmaps were generated by using package “pheatmap” in R software

Survival analysis and prognostic implications of metabolic subtypes

K–M survival plots showed that patients who belong to glycan-HCC subgroup had inferior OS when compared with patients in lipid-HCC subgroup in TCGA (HR = 1.699, 95% CI 1.126 to 2.564, p = 0.004; Fig. 3a), ICGC (HR = 3.069, 95% CI 1.465 to 6.432, p < 0.001; Fig. 3b) and CPTAC (HR = 2.514, 95% CI 1.380 to 4.581, p < 0.001; Fig. 3c) cohorts (Table 2). RMST analysis showed that the RMST for OS was shorter for glycan-HCC at 1, 3, 5 years when compared with lipid-HCC (Table 3).

Fig. 3

Clinical parameters and genomic alterations correlates of metabolic subtypes. a–c Kaplan–Meier survival analyses demonstrating poorer overall survival in glycan-HCC versus lipid-HCC in TCGA (a), ICGC (b), and CPTAC (c) cohorts. Survival plots were generated by using R package “survminer”. d Heatmap showing associations between metabolic subtypes, clinical parameters, and mutational profiles. e Genomic instability metrics (aneuploidy/CNA/HRD scores) were elevated in glycan-HCC, and tumor mutation burden showed no significant difference

Table 2 Survival Probabilities Difference Between Glycan-HCCs and Lipid-HCCs

Table 3 RMST Ratio Between Glycan-HCCs and Lipid-HCCs

Multi-omics altered between glycan and lipid subgroups

We found gene mutation status among the two subgroups showed significant difference (Fig. 3d). In the discovery database, the frequency of TP53 mutation was significantly higher in the glycan-HCC subgroup (50/122, 40.98%) compared to the lipid-HCC subgroup (58/236, 24.58%) (p = 0.001). Additionally, the mutation rate of CTNNB1 was marked lower in glycan-HCC subgroup (22/122, 18.03% vs. 72/236, 30.51% p-value = 0.011). Similarly, in the ICGC database, TP53 mutations were more prevalent in the glycan subgroup (29/59, 49.15%) than in the lipid subgroup (56/170, 32.94%) (p = 0.026). For CTNNB1, the glycan subgroup exhibited a lower mutation frequency (14/59, 23.73%) compared to the lipid subgroup (74/170, 43.53%) (p = 0.007). However, the CPTAC database did not observe the mutation frequency of TP53 (26/48, 54.17% vs. 67/111, 60.36%; p = 0.467) and CTNNB1 (7/48, 14.58% vs. 23/111, 20.72%; p = 0.364) between the glycan subgroup and lipid subgroups. We also compared four genomic instability related scores between metabolic subtypes and found that aneuploidy score, CNA fraction altered, and HRD scores were higher in glycan-HCC than those in lipid-HCC. And non-silent mutation load did not show different between two subgroups (Fig. 3e). The top 10 most mutated genes for glycan-HCC (Fig. 4a) and lipid-HCC (Fig. 4b) in TCGA database were also showed different, although global mutation types were similar between two subgroups (Fig. 4c, d). We also analyzed TP53 and CTNNB1 mutation status between two subtypes (Fig. 4e, f).

Fig. 4

Mutational landscape comparisons. Top 10 mutated genes in glycan-HCCs (a) and lipid-HCCs (b) respectively. Most frequent mutated genes in glycan-HCCs and lipid-HCCs were TP53 and CTNNB1 respectively. Transition/transversion profiles for glycan-HCCs (c) and lipid-HCCs (d). Lollipop plots depicting amino acid alterations in TP53 (e) and CTNNB1 (f) across subtypes

Transcriptomics levels alterations also analyzed by determining differentially expressed genes between glycan and lipid-subtypes were determined by using Wilcoxon test in TCGA (Fig. 5a), ICGC (Fig. 5b), CPTAC (Fig. 5c) and in-house (Fig. 5d) cohorts respectively. Then, GSEA analysis also showed altered pathways between subtypes. Interestingly, we found pathways between glycan and lipid subtypes were significant consistency in these cohorts. Pathways enriched in glycan-HCC mainly were proliferation-related pathways. For example, “cell cycle” and “ribosome” rank the top five significant pathways in the glycemic group in all four cohort (Fig. 5e–h). For pathways up-regulated in lipid subgroups, most pathways belong to metabolism-related pathways. These results indicated that metabolic subtypes have good consistency across different cohorts.

Fig. 5

Transcriptomic characterization of metabolic subtypes. Volcano plots showed differentially expressed genes between glycan-HCCs and lipid-HCCs in TCGA (a), ICGC (b), CPTAC (c) and in-house cohort (d). Gene set enrichment analysis showed pathways enriched in glycan-HCC in TCGA (e), ICGC (f), CPTAC (g) and in-house cohort (h). “Cell cycle” and “Ribosome” pathways were noticed up-regulated in glycan-HCCs in all four cohorts

TIME status showed that glycan-HCCs often exhibit higher levels of immune cell infiltration (Fig. 6a–d). However, we also noticed that immune cell exhaustion-related marker score was significantly upregulated in glycan-HCCs (Fig. 6e–h).

Fig. 6

The distribution of immune cells in different subtypes. CD8 T cells, NK cells, Monocytic lineage, Myeloid dendritic cells, Neutrophils, Cytotoxic lymphocytes, B lineage, Fibroblasts, T cells and Endothelial cells in TCGA (a), ICGC (b), CPTAC (c) and in-house cohort (d) were calculated and compared between glycan-HCCs and Lipid-HCCs. The infiltration levels of immune cells were calculated by using Microenvironment Cell Populations-counter (MCP-counter) software. Exhaustion score of T cells was higher in glycan-HCCs when compared with lipid-HCCs in TCGA (e), ICGC (f), CPTAC (g) and in-house cohort (h)

Clinically practical method for the metabolic subtyping of HCCs

Considering the TCGA cohort as the discovery cohort, the top 10 genes that were specifically up-regulation expressed in glycan HCC (PKM2, MFSD10, RCC2, TEAD2, PNMA1, LMNB2, FAM117B, TMEM51, MARCKSL1, FAM60A) and lipid HCC (SERPINC1, SLC10A1, PCK2, DAO, ACSM2A, PFKFB1, SLC27A5, GLYAT, GLYATL1, ACSM2B) (Supplementary Table S4) were eventually included for glycan-lipid polarization score factors. The glycan-lipid polarization score was defined as the arithmetic average of z-score normalized glycan-specific genes and the value of lipid-specific genes, whereby a higher score represents glycan subgroup and vice as the lipid subgroup. AUC values of the glycan-lipid polarization score in all four datasets showed high distinguishing performance (AUC > 0.9) for glycan and lipid subtypes determinations (Supplementary Figures S3A–D). Calibration curves also showed glycan-lipid polarization scores were in good agreement with the actual metabolic subtypes (Supplementary Figures S3E–H).

We analyzed the metabolic subtypes of HCC and CEUS LI-RADS patterns. The results showed that the metabolic subtypes of HCCs were significantly different among the different CEUS LI-RADS patterns; that is, the percentage of glycan-HCCs were higher in the LR-M than in the LR-5 and LR-4 categories (p < 0.05) (Table 1; Fig. 7).

Fig. 7

Contrast-enhanced ultrasound (CEUS) Liver Imaging Reporting and Data System (LI-RADS). CEUS examination obtained in a 35-year-old male patient with HCC and categorized into LI-RADS 5. a A hypoechoic lesion was detected by conventional ultrasound, and CEUS showed b arterial phase homogeneous hyperenhancement, c portal phase isoenhancement and d mild and late washout in the late phase. CEUS examination obtained in a 32-year-old female patient with HCC and categorized into LI-RADS M. e A slightly hypoechoic lesion was detected by conventional ultrasound, and CEUS showed f arterial phase homogeneous hyperenhancement followed by g early washout at 30 s and h mild and late washout in the late phase

Radiomics technology was also used for distinguishing glycan-HCCs and lipid-HCCs. Radiomics study protocol was shown in Fig. 8a. Then, 107 radiomics features were extracted in in-house (Fig. 8b) and TCGA (Fig. 8c) cohorts. A radiomics model composing 19 features were developed by using LASSO algorithm (Fig. 8d–f). Radiomics model showed moderated performance for distinguishing glycan-HCCs and lipid-HCC in TCGA cohort (AUC = 0.759, 95% CI 0.602–0.915) and in-house cohort (AUC = 0.840, 95% CI 0.756–0.924) (Fig. 8g; Table 4).

Fig. 8

Radiomics algorithm for Glycan-HCCs and Lipid-HCCs distinguishing. a Workflow of radiomics analysis, includes tumor segmentation, features extraction, model development and clinical application. Heatmaps show radiomics features profiles in in-house cohort (b) and TCGA cohort (c). Heatmaps were generated by using package “pheatmap” in R software. d A total of 19 radiomics features were finally included. e LASSO coefficient profiles of the features. f Tuning parameter (λ) selection in the LASSO model used ten-fold cross-validation via minimum criteria. g ROC curves of radiomics models in training (in-house) and test (TCGA cohort)

Table 4 Prediction performances of radiomics and AFP value for distinguishing HCC metabolic subtypes

Considering differentially expressed of serum AFP levels between glycan and lipid subgroups, the diagnostic performance of AFP was also evaluated. Serum AFP had an ROC-AUC of 0.803 (95% CI 0.773–0.817) in the in-house cohort and 0.792 (95% CI 0.776–0.807) in the CPTAC cohort (Supplementary Figure S4; Table 4).

SPP1 + macrophage is elevated in glycan-HCCs

In GSE151530 dataset, we only included HCC samples that obtained pretreatment. Then, 85 metabolic pathways enrichment score was calculated for each cell. To further determine cell communications status, we extracted 14 samples that extracted tumor cells > 0 for analysis (Fig. 9a–c). Then, enrichment scores of 85 metabolic pathways for each sample were determined by average score of tumor cells. We determined the 14 samples into glycan and lipid-HCCs based on our metabolic subtyping algorithm. Next, we compared the communication intensities strength between glycan and lipid subtypes. Results demonstrated that the connections between macrophages cells were significantly enhanced in the glycan subtype (Fig. 9d). At the global signaling pathway level, the predominant signaling patterns in the glycan subtype encompassed SPP1, FN1, LAMININ pathways (Fig. 9e, f).

Fig. 9

The communication signals differed between glycan-HCC and lipid-HCC patients. Unified manifold approximation and projection (UMAP) representation of single-cell gene expression in different patients (a), cell types (b) and metabolic subtypes (c). Umap plots were generated by using “Seurat” software. d Heatmap shows differential interaction strength among different cell types. The communication proportions (e) and levels (f) of signals differed between glycan-HCC and lipid-HCC patients. Figure were generated by using “CellChat” package in R softeare

We also delineated the ligand-receptor interactions of two subtypes, signaling changes revealed that specific increased signaling for glycan (Fig. 10a) and lipid (Fig. 10b) HCCs. Results showed SPP1 was significantly up-regulated in glycan subgroup specific macrophages (Fig. 10c) and signals emitted by macrophages in the glycan subgroup are predominantly received by CD8 + T cells, especially via the receptors CD44 (Fig. 10a, d).

Fig. 10

Identification of major signaling changes in different subtypes. Bubble plots illustrate the increased (a) and decreased (b) signaling pathways in glycan-HCC, depicting the probabilities of ligand-receptor interactions. c The gene expression of SPP1 was significantly up-regulated in macrophages associated with glycan-HCCs. d The SPP1-CD44 ligand-receptor interactions (T cells —> TAMs) were markedly up-regulated in glycan-HCCs compared to those in lipid-HCCs

Discussion

Considering liver is the major metabolic organ, it is not surprising that HCCs are characterized by metabolic alterations [22]. However, markedly heterogeneous in HCC metabolic reprogramming remain unclear, which required comprehensive analysis to broaden our knowledge to this fatal malignancy. Here, glycan-HCC and lipid-HCC subtypes were determined and multi-omics data integration revealed subtypes specific molecular alterations. Additionally, we developed surrogate methods based on gene expression, non-invasive ultrasound examination and AFP results for metabolic subtypes, these methods could reduce the cost and accelerate clinical transformation for metabolic molecular subtypes. Further analyses based on single-cell RNA-seq demonstrated that high infiltration levels of SPP1 + macrophages are the major features responsible to TIME alterations in glycan-HCCs.

We found that alterations of metabolic pathways mainly stratified patients into two distinct subgroups, included glycan-HCC and lipid-HCC. TP53 and CTNNB1 mutation are the main mutational characteristics of glycan-HCC and lipid-HCC. Previous studies have reported that TP53 mutations are involved in a number of metabolic regulatory pathways, such as the balance of glycolysis and oxidative phosphorylation, limiting the production of reactive oxygen species, and contributing to the adaptation of cells to mild metabolic stress [23, 24]. Senni et al. found that CTNNB1 mutation would have a significant impact on HCC tumor metabolism, mainly reflected in enhanced fatty acid metabolism and glutamine metabolism, and weaken the dependence of tumor cells on glycolysis [25].

Considering molecular subtypes for clinical practice still remain significant challenges, including high costs, complex technological procedures, and potential batch effects in gene expression profiling [26]. Therefore, more streamlined methods are imperative for boosting routine clinical practice. Our study provided some alternative methods for HCC’s metabolic subtyping. Some molecular biomarkers could replace the results of high-throughput sequencing to some extent. We found that SLC10A1 was upregulated in lipid-HCC, which is a key gene in bile acid co-transporter [27]. PKM2 was specific expressed in glycan-HCC, Martin SP et al. determined that elevated PKM2 is associated worse survival and treatment resistance to TACE [28]. It has been reported PKM2 is upregulated in various types of cancer and may responsibly for cancer progression and metastasis. Some previous studies attempted to determine relationships between CEUS LI-RADS category and clinicopathological features of HCC, such as tumor differentiation and Ki-67 expression [29, 30]. Furthermore, a recent study has shown that LI-RADS M HCC showed inferior clinical outcome compared with LI-RADS 3–5[31]. However, the roles of CEUS LI-RADS for molecular characteristics capture are still need explored. Metabolic reprogramming significantly influences angiogenesis and vascular distribution of HCC. Metabolites such as lactate and ketones directly alter endothelial cell metabolism, leading to disorganized vascular structures [32]. Therefore, metabolic subtypes may exhibit heterogeneity in vascular status, which could influence CEUS LI-RADS classification, necessitating further research for validation. Another non-invasive technique, radiomics, has increasingly been applied to the precise diagnosis and treatment of cancer [33, 34]. The use of radiomics has become a crucial method for capturing the molecular heterogeneity of tumors, and this approach has also demonstrated potential in the current study. We also determined serum AFP could be used as a biomarker for metabolic subtypes determination. Serum AFP is the most widely serum biomarker for HCC screening and diagnosis. A recent study also showed that worse survival of AFP-positive HCCs were responsible to some molecular characteristics, such as antigen processing, interferon-γ response and immune distortion [35]. Our study not only provided a biomarker for distinguish glycan-HCC and lipid-HCC, but also broaden molecular understanding of AFP-positive HCC.

For TIME analysis, we found that lipid-HCC exhibited lower levels of immune infiltration, primarily characterized by an “immune desert” phenotype. In contrast, glycan-HCC showed higher levels of immune cell infiltration but also showed high exhausted score of immune cells. Therefore, the immune subtypes of glycan-HCC need to be further categorized, such as exhausted and activating immune type. For lipid-HCC, efforts to enhance immune cell infiltration may improve the efficacy of immunotherapy. In glycan-HCC, immunotherapy strategies could focus on reversing the exhaustion of immune cells. By analyzing single-cell RNA-seq data of HCC, we found that SPP1 + macrophages were highly enriched in glycan-HCC compared to lipid-HCC. SPP1 + macrophages have been determined in multiple cancer types and were associated with poor survival and treatment response [36, 37].

However, some limitations in our study should be noticed. First, retrospective nature of this study limited its clinical application. Future prospective validation and applicability across different populations will provide more solid results. Second, although we have proposed multiple high-throughput alternatives for subtyping, the development of a robust multi-modal predictive subtyping approach remains an area for future research. Third, the image quality of CEUS may be affected by depth of the lesion position and obesity, which can affect the classification. The combination of CT/MR LI-RADS and CEUS LI-RADS for metabolic subtypes determination could be further explored in the future. Finally, our study lacks in vitro and in vivo mechanistic validation of the key findings. For instance, the specific regulatory mechanisms underlying the role of SPP1 + macrophages remain unexplored.

Taken together, our current work provides a comprehensive and integrated analysis of HCC metabolic reprogramming using multi-cohorts’ data. Surrogate clinically practical methods provide a bridge between high-throughput sequencing molecular data and clinical applications. Single-cell RNA-seq determined key immune cells under metabolic alterations and provided novel insights for the combination of immunotherapy and metabolism. Our study not only provides metabolic biological insights into HCC but also offers prospects for translational applications in clinical work.