
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13358-0
10.1016/j.heliyon.2024.e37327
e37327
Research Article
Development and validation of a novel immune‒metabolic-Based classifier for hepatocellular carcinoma
Zhang Wenda a1
Zhou Xinyi a1
Lin Lili b1
Lin Anqi a1
Cheng Quan cd
Liu Zaoqu ef
Luo Peng luopeng@smu.edu.cn
a⁎
Zhang Jian zhangjian@i.smu.edu.cn
a⁎⁎
a Department of Oncology, Zhujiang Hospital of Southern Medical University, Guangzhou, China
b Guangdong Provincial Key Laboratory of Malignant Tumor Epigenetics and Gene Regulation, Medical Research Center, Sun Yat-sen Memorial Hospital of Sun Yat-sen University, Guangzhou, China
c Department of Neurosurgery, Xiangya Hospital, Central South University, Changsha, Hunan, China
d National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
e State Key Laboratory of Proteomics, Beijing Proteome Research Center, National Center for Protein Sciences (Beijing), Beijing Institute of Lifeomics, Beijing, China
f Institute of Basic Medical Sciences, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
⁎ Corresponding author. luopeng@smu.edu.cn
⁎⁎ Corresponding author. zhangjian@i.smu.edu.cn
1 Joint Authors. These authors contributed equally to this work and share first authorship.

02 9 2024
15 9 2024
02 9 2024
10 17 e3732727 1 2024
4 8 2024
1 9 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
The heterogeneity of immune cells and metabolic pathways in hepatocellular carcinoma (HCC) patients has not been fully elucidated, leading to diverse clinical outcomes. Accurately distinguishing different HCC subtypes and recommending appropriate treatments is are highly important. In this study, we conducted a comprehensive analysis of 28 immune cells and 85 metabolic pathways in the TCGA-LIHC and GSE14520 datasets. Metabolism-related first principal component (MRPC1) and cytotoxic T lymphocyte (CTL) infiltration were used to assess the metabolic and immune infiltration levels of HCC patients, respectively. These two quantifiable indicators were then used to construct an immune‒metabolic classifier, which categorized HCC patients into three distinct groups. The potential biological mechanisms were explored through multiomics analysis, revealing that group S1 exhibited high metabolic activity and a high level of immune infiltration, that group S2 presented a low level of immune infiltration, and that group S3 presented low metabolic activity. This new immune‒metabolic classifier was well validated in a pancancer cohort of 9296 patients. The efficacy of multiple treatment approaches was assessed in relation to different immune‒metabolic groups, indicating that group S1 patients may benefit from immunotherapy, that group S2 patients are suitable for transcatheter arterial chemoembolization (TACE), and that group S3 patients are appropriate candidates for tyrosine kinase inhibitors. In conclusion, this immune‒metabolic classifier is anticipated to address the differences in treatment efficacy among HCC patients due to the heterogeneity of the tumor microenvironment, and to help refine the individualized treatment choices for clinical patients.

Keywords

Hepatocellular carcinoma
Tumor heterogeneity
Immune cells
Metabolism
Immunotherapy
Prognosis
==== Body
pmc1 Introduction

Hepatocellular carcinoma (HCC) accounts for approximately 90 % of liver malignancies and ranks fourth among the causes of cancer-related deaths [1]. Among the main treatments for HCC patients, transcatheter arterial chemoembolization (TACE) plays an important role, and the development of targeted therapy and immunotherapy has led to promising benefits in the treatment of HCC [2,3]. However, tumor heterogeneity is a significant characteristic of HCC [4]. Owing to the heterogeneity of immune cells and metabolic pathways in HCC patients, it is difficult to accurately distinguish patients who respond to targeted therapy or immunotherapy [5,6]. As of August 2022, the 5-year survival rate for patients with cancer of the liver and intrahepatic bile ducts was only 20.8 % [7].

Determining the appropriate stratification and effective classification of molecular subtypes can help in the development of individualized treatment strategies and improve the clinical efficacy of treatment for HCC patients [8]. Researchers have used different molecular features to define different HCC subtypes by means of CpG sites [9], protein expression profiles [10], and multiomics data [11], among other methods. However, these HCC classification models, which were developed by unsupervised or machine learning, lack interpretability, which increases the uncertainty of the models in clinical practice. Therefore, it is necessary to find suitable quantitative indicators, and establish an interpretable classifier for clinical practice.

Immune cells are involved in the immune response, immune escape, cancer infiltration and metastasis [12]. For example, activated DCs present antigens to activate T cells [13], and Treg cells play a prominent immunosuppressive role in blocking the anticancer immune response [14]. Cytotoxic T lymphocyte (CTL) are capable of secreting cytokines and killing tumor cells, which is a key role of defense in antitumor immunity [15]. Moreover, metabolic reprogramming is one of the hallmarks of human cancer [16]. Immune cells and metabolic pathways play diverse and important roles in the development of HCC. Therefore, researchers are attempting to establish metabolic or immune types through different metabolic pathways or immune cell compositions in HCC patients [17,18]. However, just as the immune cell-activating energy supply switches from oxidative phosphorylation to aerobic glycolysis, glutamine deficiency inhibits macrophage and T-cell functions [19,20], immune cells and metabolic pathways do not exist independently of each other. The degree of immune cell infiltration and activation of metabolic pathways are important components of heterogeneity in the tumor microenvironment (TME) of HCC patients [21,22]. Therefore, combining immune cell infiltration and metabolic pathway activity to develop an “immune‒metabolic classifier” will be an important direction for establishing novel molecular subtypes of HCC. However, how to accurately quantifying and combining immune cells and metabolic pathways is one of the current difficulties.

In this study, we aimed to identify quantifiable indicators of immune cells and metabolic pathways, with the goal of establishing a novel immune‒metabolic classifier for predicting molecular subtypes in HCC patients. This classifier was developed through a systematic analysis of 85 metabolic pathways and 28 immune cell types [23,24]. The validity of the classifier will be assessed via multiomics analysis to explore the underlying biological mechanisms. Furthermore, we investigated the potential beneficial treatment options for HCC patients across different immune‒metabolic subtypes, aiming to provide guidance for personalized treatment strategies. In summary, this immune‒metabolic classifier is anticipated to address the differences in treatment efficacy among HCC patients due to the heterogeneity of the TME and to help refine the individualized treatment choices for clinical patients.

2 Materials and methods

This study adhered to the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD) guidelines. The primary endpoint was overall survival (OS). The study protocol was approved by the ethics committee of Zhujiang Hospital. The data used in this study were obtained from publicly available datasets, and the requirement for informed consent was waived.

2.1 Patients and datasets

Bulk RNA-sequencing (RNA-seq) data from 10,779 patients including 834 HCC patients, 9296 pancancer patients and 649 patients with a history of immunotherapy, were included in this study. The patients included in the study had detailed survival information or treatment information. In addition, single-cell RNA sequencing (scRNA-seq) data from ten HCC tissues were extracted from the GSE149614 dataset. The detail HCC patient information is provided in Appendix Table S1. GEO datasets were obtained from GEO (http://www.ncbi.nlm.nih.gov/geo/). TCGA datasets were obtained from the UCSC Xena browser (GDC hub: https://gdc.xenahubs.net). IMvigor210 and PRJEB23709 were obtained from published literature [25,26]. All the data were collected from March 25 to May 27, 2022. Basic patient information is provided in Appendix Table S2.

2.2 Data preprocessing

For the bulk RNA-seq data, the FPKM values were transformed into TPM values. For the microarray data, the normalizeBetweenArrays function in the limma R package was used to normalize the data.

For the scRNA-seq data, the Seurat R package (v 4.0.0) [27] was used for quality control of individual cells, data normalization, unsupervised clustering, and visualization. The selection criteria were described in previous studies [28] and low-quality cells (<200 genes/cell or >8000 genes/cell and >20 % of mitochondrial genes) were removed. The top 2000 variable features were identified for each sample, and integration was achieved via canonical correlation analysis (CCA) algorithms. Principal component analysis (PCA) was used to reduce the dimensionality, and the top 50 principal components were selected for UMAP.

For the methylation data, the ChAMP R package (v2.24.0) [29] was used for quality control and analysis. The champ.filter function was used to filter out low-quality probes and to remove samples with NA values. The function of CHAM.DMP was used to detect differentially methylated genes with the criteria of |log2FC|>0.4 and adjusted P < 0.05 [30].

For the somatic mutation data, the oncoplot function in the maftools R package (v2.20.0) [31] was used to identify and visualize the frequently mutated genes.

2.3 Clustering for immune cells and metabolic pathways

The GSVA R package (v1.42.0) [32] was utilized for the quantification of immune cells and metabolic pathways in the TCGA-LIHC and GSE14520 datasets. The ConsensusClusterPlus R package (v1.58.0) [33] was used to identify molecular subtypes of HCC on the basis of the activity of immune cells and metabolic pathways. At each iteration 80 % of the samples were divided into subsamples and each subsample was divided into up to k (maximum K = 6) groups basis of the on the Euclidean distance via the k-means algorithm. This process was repeated 1000 times. The optimal number of clusters was determined by cumulative distribution function (CDF) curves, proportion of ambiguous clustering (PAC) scores and Nbclust testing, which should have the lowest PAC score and be supported by the highest number of criteria. The same clustering analysis was performed on the TCGA-LIHC and GSE14520 datasets to ensure the stability of the clustering results.

2.4 scRNA-seq data analysis

The FindCluster function (resolution = 0.5) was employed to identify cell clusters. Cell-specific markers were used to identify cell types, and these markers were obtained from a previous study [21] and the CellMarker database [34]. The AUCell package (v1.16.0) [35] in R was used to estimate the gene set activity of each cell, and these gene sets were generated with the Molecular Signatures Database (MSigDB) of the Broad Institute [36]. The CellChat R package (v1.6.1) [37] was used to identify the patterns of cell‒cell communication.

2.5 Identifying the quantifiable indicators of metabolic pathways

The DESeq2 R package [38] was used to identify (DEGs) with the criteria of |log2foldchange| > 2 and P < 0.05. Metabolism-related differential genes (MRDEGs) were defined as the intersection of DEGs and metabolic genes. metabolism-related first principal component (MRPC1) was defined as the first principal component of the MRDEGs, and was used to assess the activity of metabolic pathways.

2.6 Internal mechanism exploration

The metabolic pathways and immune-related signatures were used to assess the immune‒metabolic landscape in the different groups; the related signatures are shown in Appendix Table S3. The clusterProfiler R package (v4.0.5) [39] was utilized for enrichment analysis and visualization, including Kyoto Encyclopedia of Genes Genomes (KEGG), Gene Ontology (GO) and and gene set enrichment analysis (GSEA). Pathways were considered significantly different when the P value and false discovery rate were less than 0.05. Coexpressed module eigengenes were identified via the WGCNA R package [40]. The random forest algorithm was used to rank the importance of signatures and was implemented through the randomForest R package (v4.7-1.1). The level of CTL infiltration was calculated on the basis of the average expression of PRF1 and GZMA.

2.7 Construction of a novel immune‒metabolic classifier

A novel immune‒metabolic classifier was constructed with the optimal cutoff values for MRPC1 and CTL infiltration in the TCGA-LIHC cohort, which were used to predict the three immune‒metabolic groups, designated as S1, S2 and S3. The immune‒metabolic classifier first distinguished patients into the S3 group and the combined S1/S2 groups on the basis of the optimal cutoff value of MRPC1, and then the S1/S2 group was distinguished into the independent S1 and S2 groups on the basis of the optimal cutoff value of CTL infiltration. The GSE14520 dataset was used as a parallel verification cohort to ensure the stability of the analysis. Three other HCC datasets, IMvigor210 and a pancancer cohort containing 32 tumor types were also used to test the stability of the classifier. To avoid the potential batch effects, the immune‒metabolic classifiers was used individually for each tumor type in the pancancer cohort, and their prediction results were then combined. The surv_cutpoint function of the survminer R package (v0.4.9) was used to establish the optimal cutoff value for continuous variables (if survival data were missing from the dataset, the median value was taken as the cutoff value). To assess the classifier prediction performance, the confusionMatrix function of the caret R package (v6.0-94) was applied to evaluate its accuracy, sensitivity, recall, and F1 score for each group, which was subsequently weighted with the percentage of patients in each group to calculate the weighted sensitivity, weighted recall, and weighted F1 score.

2.8 Clinical treatment evaluation

The pRRophetic R package (v0.5) [41] is based on the Genomics of Drug Sensitivity in Cancer (GDSC) database [42] and was used to evaluate the half-maximal inhibitory concentration (IC50) of each group of therapeutic agents to screen for potential therapeutic agents.

2.9 Constructing a novel immunotherapy predictor

Weighted gene coexpression network analysis (WGCNA) and gene enrichment analysis were employed to identify of immune-related module eigengenes in the TCGA-LIHC cohort. The S1 signature represents the intersection of DEGs with immune-related module eigengenes and was further applied to construct the S1 score via the single sample gene set enrichment analysis (ssGSEA) algorithm. TCGA-LIHC, GSE14520, and a pancancer cohort containing 32 tumor types were used to assess the prognostic performance of the S1 score. Seven immunotherapy cohorts were used to assess the benefit of immunotherapy in different S1 score groups.

2.10 Statistical analysis

All the statistical analyses and visualizations were performed via R (version 4.1.2). Differences in OS and progression-free survival (PFS) between different subgroups were evaluated via survival curves. Hazard ratios (HRs) and 95 % confidence intervals (CIs) were assessed via univariate Cox proportional hazards regression models. The Pearson coefficient was used to measure the correlation between two continuous variables, the chi-square test was used to compare two categorical variables, and the Shapiro‒Wilk normality test was used to check the normality of the distribution of variables. For normally distributed variables, unpaired Student's t tests were used to compare two groups; one-way analysis of variance (ANOVA) was used for more than two groups. For nonnormally distributed variables, the Wilcoxon rank-sum test was used to compare two groups; the Kruskal‒Wallis test was used for more than two groups. Differential expression analysis was performed via the DESeq2 package in R to identify DEGs. A multivariate Cox regression model was applied to identify independent prognostic factors. Finally, receiver operating characteristic (ROC) curves were constructed to evaluate the sensitivity and specificity of the selected features.

3 Results

3.1 Identification of two immune‒metabolic clusters with different clinical outcomes

To explore the differences in the prognosis of HCC patients due to the immune cell and metabolic pathway heterogeneity, 28 infiltrating immune cells and 85 metabolic pathway activities were assessed in the TCGA-LIHC cohort. The HCC patients were divided into two clusters (C1 and C2) via consensus clustering (Fig. 1A). Compared with cluster C2, cluster C1 had a significantly greater survival benefit than cluster C2 (OS: HR = 1.70; 95 % CI, 1.17–2.47; P = 0.0057; PFS: HR = 1.73; 95 % CI, 1.18–2.52; P = 0.0046) (Fig. 1B). These findings were verified in HCC patients from the GSE14520 cohort (OS: HR = 2.28; 95 % CI, 1.49–3.49; P = 0.0001; PFS: HR = 1.62; 95 % CI, 1.13–2.33; P = 0.0084) (Fig. 1C).Fig. 1 Identification of two immune‒metabolic clusters of HCC with different clinical outcomes. (A) The proportion of ambiguous clustering (PAC) score, a low PAC value indicates a flat middle segment, allowing conjecture of the optimal k (k = 2) by the lowest PAC value. Number of the criteria for the recommended number of clusters were identified via the NbClust package. (B, C) Kaplan‒Meier curves of overall survival and progression-free survival based on the C1 and C2 clusters in the TCGA-LIHC and GSE14520 cohorts, respectively. (D) The immune and metabolic patterns of the two clusters. (E) The random forest method was used to rank the importance of 16 hallmark pathways. The more important the hallmark pathway is, the greater the mean decrease in the Gini coefficient is. (F) Overlap of DEGs and metabolic genes. (G) ROC curves showing the ability of MRPC1 to predict two clusters of HCC. DEGs, differentially expressed genes; MRPC1, metabolism-related first principal component; ROC, receiver operating characteristic.

Fig. 1

The differences between the two clusters were related mainly to metabolic pathway activities (Fig. 1D– Appendix Fig. S2A). GSEA revealed 16 significant pathways, the bile acid pathway being the most important metabolic pathway that with the highest MeanDecreaseGini value (Fig. 2H; Appendix Fig. S2B, C). High bile acid pathway activity was associated with improved OS (Appendix Fig. S2D), and cluster C2 expressed a low level of SIRT5, which degrades primary bile acids (Appendix Fig. S2E). These findings support the classification of HCC patients into C1 and C2.Fig. 2 Comparison of cell‒cell communication between the C1 and C2 clusters (A) Nine major cell types in HCC tissues. (B) Proportion of diverse cell types in the two clusters. (C) Distribution of MRPC1 and eight metabolic pathways (D) Expression of representative markers across diverse cell types. (E) The number and strength of ligand–receptor pairs between the C1 and C2 clusters. (F) Differential expression analysis at the pathway level between the C1 and C2 clusters. (G) Activity of 25 antigen presentation-related pathways in the two clusters. (H) All significant signaling pathways are ranked according to the difference in overall information flow between the C1 and C2 clusters. (I) Communication probability of ligand–receptor pairs of DCs and other cell types. DC, Dendritic cell; HSC, Hematopoietic stem cell; MRPC1, Metabolism-related first principal component.

Fig. 2

To investigate the effects of different metabolic patterns on the interactions of various cell types in the TME, differential expression analysis was performed for the two clusters in the TCGA-LIHC and GSE14520 cohorts. A total of 134 MRDEGs were obtained and enriched in nine metabolic pathways. Among these pathways, eight of which were highly activated in cluster C1 (Fig. 1F–Appendix Fig. S3A, B). These MRDEGs were subsequently used for principal component analysis and the first principal component was taken as MRPC1. MRPC1 was subsequently used as a quantifiable index to distinguish the C1 and C2 clusters from the TCGA-LIHC and GSE14520 cohort, with ROC values of 0.93 and 0.82, respectively (Fig. 1G).

Furthermore, the changes in the TME of HCC patients at different metabolic levels were investigated by using the scRNA-seq data of HCC patients. Nine major cell types were defined by cell-specific marker expression (Fig. 2A, B, D; Appendix Fig. S3C), and 10 HCC single-cell samples were categorized into high- and low-MRPC1 groups according to the median value of MRPC1 (Fig. 2C). The high-MRPC1 group was considered the C1 cluster, whereas the low-MRPC1 group was considered the C2 cluster. The communication intensity and number of ligand‒receptor pairs between various cell types were greater in the C1 clusters than in the C2 clusters (Fig. 2E). Pathways, including TNF, MHC-I, and MHC-II, exhibited a prominent changes in information flow in the C1 cluster, whereas the NOTCH, EPHA, and EPHB pathways exhibited prominent changes in information flow in the C2 cluster (Fig. 2H). Metabolic pathways play important roles in dendritic cell (DC) differentiation, function, and proliferation, and the communication mode between DCs and other cell types differed between the C1 and C2 clusters (Fig. 2I). The pathway related to antigen presentation was more highly activated in the C1 cluster (Fig. 2F and G).

3.2 Identification of three immune‒metabolic groups with different molecular characteristics

Moreover, there were significant differences in immune cell infiltration between the two clusters. Notably, a large fluctuation in immune cell infiltration was observed in cluster C1 (Fig. 3A), and CTL infiltration was used as a quantifiable index to represent 28 immune cells because a significant positive correlation between CTL infiltration and infiltration of 28 immune cell types was observed (Fig. 3B). Cluster C1 was further divided into two subgroups according to CTL infiltration levels. The high-CTL infiltration subgroup had significantly better survival benefits than the low-CTL infiltration subgroup did (TCGA-LIHC: HR = 2.07, 95 % CI: 1.23–3.38, P = 0.0036; GSE14520: HR = 1.93, 95 % CI: 1.04–3.56, P = 0.0368) (Fig. 3C). The high-CTL infiltration subgroup was subsequently referred to as group S1 and the low-CTL infiltration subgroup was group S2; C2 was considered group S3. Among the three groups, group S1 had the longest OS (Fig. 3D).Fig. 3 The C1 cluster was further classified into two groups via CTL infiltration (A) Infiltration of 28 immune cells across two clusters. (B) Correlations between CTL and 28 immune cells. (C) Kaplan‒Meier curves of OS based on the high- and low-CTL infiltration groups in the TCGA-LIHC and GSE14520 cohorts. (D) Kaplan‒Meier curves of OS based on the three immune‒metabolic groups in the TCGA-LIHC and GSE14520 cohorts. CTL, cytotoxic T lymphocyte. OS, overall survival; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Fig. 3

To further explore the underlying molecular mechanisms involved, differences in the TME among the three groups were evaluated. Different immune‒metabolic patterns among the three groups was depicted, which revealed that group S1 exhibited high metabolic activity and a high level of immune infiltration, that group S2 presented a low level of immune infiltration, and that group S3 presented low metabolic activity (Fig. 4A and B). We further observed that tumor tissues were less methylated than normal tissues were. Among three groups, group S3 was the most highly methylated and had the highest percentage of hypermethylated genes (Fig. 4C–E). The hypermethylated genes in group S3 were highly enriched in multiple metabolism-related pathways (Fig. 4D). Somatic mutation analysis revealed that group S3 had the highest percentage of TP53 mutations and that group S2 had the highest percentage of CTNNB1 mutations (Fig. 4F). In addition, different module eigengenes were found in the three immune‒metabolic groups and presented different functions. WGCNA analysis revealed that the yellow module, which correlates positively with groups S1 and S2, primarily enriches immune-related pathways; the blue module, correlating with groups S1 and S3, is mainly enriched in metabolism-related pathways; while the brown and turquoise modules, both correlated with group S3, are primarily enriched in proliferation- and oncogene-related pathways, respectively (Fig. 5A and B, Appendix Fig. S4 A). Moreover, Group S1 had greater interferon pathway activity and higher chemokine and chemokine receptor expression than did group S2 (Fig. 5B–D, Appendix Fig. S4B-D). Taken together, these three immune‒metabolic groups have different survival benefits and molecular characteristics.Fig. 4 Landscape of immune‒metabolic patterns, methylation alterations and somatic mutations (A, B) Immune and metabolic patterns of the three immune‒metabolic groups. (C) Levels of promoter methylation in the three immune‒metabolic groups and adjacent normal samples in the TCGA-LIHC cohort. (D) Pathway enrichment analysis of three immune‒metabolic groups of hypermethylated genes. (E) The proportion of differentially methylated CpG sites was determined using adjacent normal samples as a reference. (F) Mutation frequency of the top 20 FMGs among the three immune‒metabolic groups. FMGs; FMGs, Frequently mutated genes. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Fig. 4

Fig. 5 Identification of three immune‒metabolic groups with different molecular characteristics (A) Association of three immune‒metabolic groups, PFS, DFS and OS with MEs in the TCGA-LIHC cohort. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of four MEs in the TCGA-LIHC cohort. (C) Distribution of the six hallmark pathways between the S1 and S2 groups. (D) Gene set enrichment analysis showing DEGs between S1 and S2 groups evaluated in the context of hallmark gene sets. (E) Distribution of chemokines and chemokine receptors between the S1 and S2 groups. DFS, disease-free survival; MEs, module eigengenes; OS,overall survival; PFS, progression-free survival.

Fig. 5

3.3 An immune‒metabolic-based classifier was established and validated in multiple cancer patients

Metabolic pathways are closely related to immune cells in the TME (Fig. 6A). A classifier was generated on the basis of the cutoff values for MRPC1 and CTL infiltration, and the classification efficacy of the classifier was evaluated (TCGA-LIHC: accuracy: 0.90, weighted F1: 0.90, weighted sensitivity: 0.93, weighted recall: 0.90; GSE14520: accuracy: 0.75, weighted F1: 0.75, weighted sensitivity: 0.83, weighted recall: 0.75) (Fig. 6B–E). The patients in group S1 had significantly longer OS than those in the other two groups did (Fig. 6C). Moreover, HCC patients were divided into different clinical subgroups according to age, sex, AFP, and stage. Survival analysis revealed that patients in group S1 tended to have the best survival prognosis among the different clinical subgroups, suggesting that the classifier is partially independent of existing clinical characteristics (Fig. 6F). The classifier was applied to the pancancer cohort. Among the tumor types, group S1 patients had the highest OS, and the proportions of patients in the three immune‒metabolic groups were highly consistent (Fig. 6D–G; Appendix Fig. S5).Fig. 6 An immune‒metabolic classifier was established and validated in multiple cancer patients (A) Correlations between CTL infiltration, MRPC1 expression, the activity of eight metabolic pathways and the infiltration of 28 immune cells in the TCGA-LIHC cohort. (B) CTL and MRPC1 distributions and their expression among the three immune‒metabolic groups in the TCGA-LIHC cohort. (C, D) Kaplan‒Meier curves of OS in the three immune‒metabolic groups. (E) Performance of the immune‒metabolic classifier in predicting three immune‒metabolic groups. (F) Patients in group S1 tended to have the best survival prognosis among the different clinical subgroups. S1 was used as a reference. (G) Patients in group S1 tended to have the best survival prognosis among patients with various tumor types. S1 was used as a reference. CTL, cytotoxic T lymphocyte; MRPC1, metabolism-related first principal component; OS, overall survival.

Fig. 6

3.4 HCC patients in the three immune‒metabolic groups are suitable for different therapeutic interventions

To determine the appropriate treatment strategy for the three immune‒metabolic groups, we evaluated the survival benefit and response rate across treatments. compared with the other two groups, group S1 patients had the best survival benefits and highest response rates to immunotherapy(Fig. 7A). These findings were validated in patients treated with the PD-1 inhibitor from IMvigor210 (Fig. 7B). Submap analysis also revealed that group S1 patients benefitted more from immunotherapy (Fig. 7C). Patients in groups S1 and S2 benefitted from TACE, but patients in group S3 were more likely to benefit from sorafenib (Fig. 7D). Patients in group S3 were more sensitive to sorafenib, nilotinib, and bosutinib according to the assessment of drug sensitivity by the IC50 (Fig. 7E).Fig. 7 HCC patients in the three immune‒metabolic groups are suitable for different therapeutic interventions (A–B) Kaplan‒Meier curves of OS based on the three immune‒metabolic groups. Objective response rate stratified by three immune‒metabolic groups. (C) Group S1 was the most sensitive to anti-PD1 therapy (Bonferroni, P = 0.011). (D) Heatmap showing the expression of MRPC1 and CTL infiltration. Bar plot showing the treatment response ratio among the three immune‒metabolic groups. (E) Distribution of half-maximal inhibitory concentration (IC50) values among the three immune‒metabolic groups. CR, complete response; CTL, cytotoxic T lymphocyte; CTLA4, cytotoxic T lymphocyte-associated protein 4; MRPC1, metabolism-related first principal component; OS, overall survival; PD, progressive disease; PR, partial response; SD, stable disease; NR, no response; PD1, programmed cell death protein 1; R, response. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Fig. 7

3.5 The S1 score can serve as a prognostic biomarker and predict immunotherapy efficacy

To further accurately distinguish which patients are suitable for immunotherapy, a prognostic model, named the S1 score, based on thirty S1 signature genes was established (Fig. 8A), and patients with HCC were divided into high-S1 score and low-S1 score groups. The high-S1 score group had longer OS than the low-S1 score group did (TCGA-LIHC: HR = 2.00, 95 % CI: 1.40–2.85, P = 0.0001; GSE14520: HR = 1.97, 95 % CI: 1.27–3.05, P = 0.0025) (Fig. 8B and C). In addition, multivariate Cox regression analysis revealed that the S1 score was an independent prognostic factor (Fig. 8D), and the predictive efficacy of the S1 score was well validated across multiple tumors in the pancancer cohort (Fig. 8E).Fig. 8 The S1 score is a prognostic biomarker and predicts immunotherapy efficacy (A) Distribution of S1 signature gene among the three immune‒metabolic groups. (B, C, F-I) Kaplan‒Meier curves of OS based on high- and low-S1 scrore groups. (D) Multivariate Cox regression analysis of OS in the TCGA-LIHC cohort. (E) Prognostic performance of the S1 score across cancers. OS, overall survival.

Fig. 8

Furthermore, the S1 score performed well in predicting the efficacy of immunotherapy. The high-S1 score group had significantly longer OS than the low-S1 score group did (PRJEB23709: HR = 3.82, 95 % CI: 1.96–7.46, P = 0.0001; IMvigor210: HR = 1.43, 95 % CI: 1.09–1.88, P = 0.0102; GSE91061-Pre: HR = 3.03, 95 % CI: 1.23–7.43, P = 0.0153; GSE91061-On: HR = 2.42, 95 % CI: 1.15–5.10, P = 0.0206; GSE140901: HR = 1.66, 95 % CI: 0.65–4.21, P = 0.2910; GSE106128: HR = 5.02, 95 % CI: 1.20–21.02, P = 0.0273) (Fig. 8F–I; Appendix Fig. S6A, B). The prediction efficiency of the S1 score for immunotherapy was also assessed (area under the curve (AUC): IMVigor210, 0.54; GSE35640, 0.72; GSE140901, 0.58; GSE91061-Pre, 0.65; GSE91061-On, 0.76; PRJEB23709, 0.81; GSE106128, 0.63) (Appendix Fig. S7C). Taken together, the results emphasize that patients in the high-S1 score group had better clinical outcomes and higher immunotherapy response rates for multiple types of tumors.

4 Discussion

Despite the excellent performance of targeted therapy [43,44] and immunotherapy in HCC [45], the heterogeneity of immune cells and metabolic pathways in HCC has largely been unexplored and has led to different clinical outcomes. In this study, the quantifiable indicators MRPC1 and CTL were identified, which used to develop a new immune‒metabolic classifier. This new immune‒metabolic classifier divides HCC patients into three groups and was well validated in the pancancer cohort. The associations of different immune‒metabolic groups with the efficacy of multiple treatment approaches were assessed, group S1 patients may benefit from immunotherapy, group S2 patients are suitable for TACE, and group S3 patients are appropriate candidates for tyrosine kinase inhibitors. Furthermore, the S1 score, a prognostic model for predicting immunotherapy efficacy in pancancer patients, was established on the basis of thirty signature genes of group S1. These results could offer insights into the development of individualized treatment strategies for HCC patients and provide valuable references for the selection of clinical therapeutic strategies.

The metabolic microenvironment is complicatedly linked to immune cell function, proliferation, and differentiation [46,47]. Tumor cells and immune cells strongly compete for metabolism [46]. On the one hand, with highly active metabolic pathways, tumor cells can easily develop deficiencies in oxygen and critical nutrients and generate lactic acid. Lactic acid accumulation leads to acidosis, thus hindering immune cell function, proliferation and differentiation [48,49]. On the other hand, the inadequate supply of energy and critical nutrients to immune cells may also be an important cause of immune dysfunction and immune escape in cancer patients. For example, the function, maturation and differentiation of immune cells in DCs also require the involvement of metabolic pathways such as glycolysis, the tricarboxylic acid cycle, and corresponding metabolic intermediates [48]. When the corresponding energy supply and nutrients are lacking, DC antigen-presenting ability is reduced, limiting T-cell function. Owing to the important role of metabolic pathways in TME, this study defines a quantifiable indicator, MRPC1, to accurately measure its status, which is the first key result of this study.

As we mentioned above, the energy and nutrient supplies of immune cells are essential factors in the response of immune cells to immunotherapy [50]. In this study, we found that the survival of patients in group S1 improved the most, which was largely due to its highly activated energy metabolism pathway, which can provide sufficient energy and nutrients for immune cells. Although patients in group S3 exhibited a greater degree of immune cell infiltration than patients in group S2 did, the energy supply in group S3 was inadequate, which may be one of the reasons why these patients had worse OS than did those in group S2. The multiomics data suggest that the hypermethylated genes in group S3 are enriched in multiple metabolism-related pathways, which may contribute to the low metabolic activity of these patients. In addition, it has been reported that reduced SIRT5 expression can cause the accumulation of primary bile acid, which may further promote M2 macrophage polarization and suppress the activity of immune cells [51]. Consistent with these findings, the metabolic microenvironment of group S3 may have inhibited immune cell function because patients in group S3 presented lower SIRT5 expression than did those in the other two groups. The different immune‒metabolic statuses presented by these three groups are also the reason for the transformation of HCC patients from two clusters to three subtypes in this study.

Patients given the same treatment may exhibit substantial differences in efficacy because of the different molecular characteristics of different HCC subtypes [52]. TACE is a first-line topical therapy for patients with intermediate or advanced HCC, and is increasingly being administered in combination with other therapies, including sorafenib and immunotherapy [53,54]. Our study revealed that patients in groups S1 and S2 were likely to benefit from TACE therapy. Patients in group S1 also benefit from immunotherapy. Patients with “immune hot” tumors are more likely to benefit from immunotherapy [55]. Consistent with this, patients in group S1 exhibited an “immune hot” tumor subtype with abundant immune cell infiltration and sufficient energy and nutrients. Patients in group S1 also had a strong antigen-presenting capacity to support immune cell function. PDL1 and CTLA4 can inhibit the function of T cells, the main participants in antitumor activity, and are targets of immunotherapy [56]. Similarly, we found that patients in group S1 had a high degree of immune cell infiltration and high immune checkpoint expression, and the characteristic genes of group S1 patients were enriched in several immune-related pathways. Therefore, we further established the S1 score on the basis of 30 characteristic genes in the group S1, and the patients in the high-S1 score group had better clinical outcomes and higher immunotherapy response rates.

Additionally, patients in group S2 had the highest percentage of CTNNB1 mutations, which is consistent with the finding that CTNNB1 mutations predispose patients toward immunotherapy resistance [57]. Patients in group S3 had the highest percentage of TP53 mutations, which is consistent with the finding that TP53 mutations are associated with immune escape and poor prognosis [58]. These features might be components of the potential mechanisms for the poor benefits of immunotherapy in patients in groups S2 and S3. The signature genes of patients in group S3 were enriched in multiple oncogenic pathways, which suggested that these patients were more likely to benefit from treatment with oncogenic pathway inhibitors. In this study, we found that patients in group S3 were sensitive to sorafenib, and these findings were confirmed via the PRRophetic algorithm.

Despite the heterogeneity of cancers in terms of multiomics data, such as transcriptomics, proteomics and genomics data, certain important immune cells and metabolic pathways play the same role in different types of cancer. For example, CTLs can play an essential role in killing tumor cells in different types of cancer, and the bile acid metabolic pathway is key as a signaling molecule in regulating metabolic and immune dynamic homeostasis in various organs throughout the body. Immune cells and metabolic pathways have crucial functions in different types of tumors. Hence, the new immune‒metabolic classifier and S1 score could be used for many other types of tumors.

This study has several limitations. First, HCC datasets included in this study lack sufficient clinical information, which may cause bias in patient selection and lack of comprehensive validation for the classifier. Senond, in vitro and in vivo experiments are needed to reveal the molecular basis of the immune‒metabolic classifier and the effectiveness of putative drugs for the treatment of HCC. Third, many targeted therapy drugs, including lenvatinib [59] and regorafenib [60], have been applied in HCC treatment and have achieved good efficacy. Rigafenib has also been reported to exert immunomodulatory effects and to synergize withimmune checkpoint inhibitors [61]. Although we suggest that patients in group S3 are likely to benefit from treatment with oncogenic pathway inhibitors, it is still unknown whether group S3 patients can benefit from those drugs or not because of the data accessibility issues. Additionally, these two predictive models still need to be validated with large samples, and multicenter datasets and a validation dataset is being collected.

5 Conclusion

This study established a new immune‒metabolic classifier and analyzed the differences in clinical outcomes, metabolic pathways, immune landscapes, and treatment responses among patients in three immune‒metabolic groups. The immune‒metabolic classifier has precise and robust efficacy prediction and can provide individualized treatment strategies for HCC patients, including TACE, targeted therapy and immunotherapy. The S1 score can serve as a prognostic biomarker for immunotherapy in pancancer patients. These two predictive models may serve as a reference for the individualized treatment choices for clinical patients, and are worth recommending.

Data availability

Data to support the findings of this study is available on reasonable request from the corresponding author. The datasets included in this study were obtain from the UCSC Xena browser (GDC hub: https://gdc.xenahubs.net) and GEO (http://www.ncbi.nlm.nih.gov/geo/).

Funding

None.

Ethics approval

None.

CRediT authorship contribution statement

Wenda Zhang: Writing – review & editing, Writing – original draft, Visualization, Software, Formal analysis, Conceptualization. Xinyi Zhou: Writing – review & editing, Writing – original draft, Software, Formal analysis, Conceptualization. Lili Lin: Writing – review & editing, Writing – original draft, Visualization, Software, Formal analysis. Anqi Lin: Writing – original draft, Visualization, Formal analysis. Quan Cheng: Visualization, Software. Zaoqu Liu: Visualization, Software. Peng Luo: Writing – review & editing, Writing – original draft, Software, Resources, Project administration, Investigation, Conceptualization. Jian Zhang: Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Methodology, Investigation, Data curation, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following are the Supplementary data to this article:Multimedia component 1

Multimedia component 1

Multimedia component 2

Multimedia component 2

Multimedia component 3

Multimedia component 3

Multimedia component 4

Multimedia component 4

Acknowledgments

None.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e37327.
==== Refs
References

1 Donne R. Lujambio A. The liver cancer immune microenvironment: therapeutic implications for hepatocellular carcinoma Hepatology 77 2023 1773 1796 35989535
2 Reig M. Forner A. Rimola J. BCLC strategy for prognosis prediction and treatment recommendation: the 2022 update J. Hepatol. 76 2022 681 693 34801630
3 Sangro B. Sarobe P. Hervás-Stubbs S. Melero I. Advances in immunotherapy for hepatocellular carcinoma Nat. Rev. Gastroenterol. Hepatol. 18 2021 525 543 33850328
4 Wang C. Qin S. Pan W. mRNAsi-related genes can effectively distinguish hepatocellular carcinoma into new molecular subtypes Comput. Struct. Biotechnol. J. 20 2022 2928 2941 35765647
5 Fang H. Guo Z. Chen J. Combination of epigenetic regulation with gene therapy-mediated immune checkpoint blockade induces anti-tumour effects and immune response in vivo Nat. Commun. 12 2021 6742 34795289
6 Wei C. Zhu M. Zhang P. PKCα/ZFP64/CSF1 axis resets the tumor microenvironment and fuels anti-PD1 resistance in hepatocellular carcinoma J. Hepatol. 77 2022 163 176 35219791
7 Devarbhavi H. Asrani S.K. Arab J.P. Global burden of liver disease: 2023 update J. Hepatol. 79 2023 516 537 36990226
8 Rebouissou S. Nault J. Advances in molecular classification and precision oncology in hepatocellular carcinoma J. Hepatol. 72 2020 215 229 31954487
9 Cheng J. Wei D. Ji Y. Integrative analysis of DNA methylation and gene expression reveals hepatocellular carcinoma-specific diagnostic biomarkers Genome Med. 10 2018 42 29848370
10 Gao Q. Zhu H. Dong L. Integrated proteogenomic characterization of HBV-related hepatocellular carcinoma Cell 179 2019 561 577.e22 31585088
11 Cancer Genome Atlas Research Network Cancer genome atlas research network, comprehensive and integrative genomic characterization of hepatocellular carcinoma Cell 169 2017 1327 1341.e23 28622513
12 Lu Y. Yang A. Quan C. A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma Nat. Commun. 13 2022 4594 35933472
13 Wculek S.K. Cueto F.J. Mujal A.M. Dendritic cells in cancer immunology and immunotherapy Nat. Rev. Immunol. 20 2020 7 24 31467405
14 Granito A. Muratori L. Lalanne C. Hepatocellular carcinoma in viral and autoimmune liver diseases: role of CD4+ CD25+ Foxp3+ regulatory T cells in the immune microenvironment, World J. Gastroenterol. 27 2021 2994 3009
15 Weber E.W. Maus M.V. Mackall C.L. The emerging landscape of immune cell therapies Cell 181 2020 46 62 32243795
16 Gong Y. Ji P. Yang Y. Metabolic-pathway-based subtyping of triple-negative breast cancer reveals potential therapeutic targets Cell Metabol. 33 2021 51 64.e9
17 Wei X. Michelakos T. He Q. Association of tumor cell metabolic subtype and immune response with the clinical course of hepatocellular carcinoma Oncol. 28 2023 e1031 e1042
18 Zhai W. Lai H. Kaya N.A. Dynamic phenotypic heterogeneity and the evolution of multiple RNA subtypes in hepatocellular carcinoma: the PLANET study Natl. Sci. Rev. 9 2022 nwab192
19 Riera-Domingo C. Audigé A. Granja S. Cheng W. Immunity, Hypoxia, and Metabolism-the Ménage à Trois of Cancer: Implications for Immunotherapy Physiol. Rev. 100 2020 1 102 31414610
20 Hoogerland J.A. Staels B. Dombrowicz D. immune‒metabolic interactions in homeostasis and the progression to NASH Trends Endocrinol. Metab. 33 2022 690 709 35961913
21 Chen D. Zhang Y. Wang W. Identification and characterization of robust hepatocellular carcinoma prognostic subtypes based on an integrative metabolite-protein interaction network Adv. Sci. 8 2021 e2100311
22 Lei X. Lei Y. Li J.-K. Immune cells within the tumor microenvironment: biological functions and roles in cancer immunotherapy Cancer Lett. 470 2020 126 133 31730903
23 Xiao Z. Dai Z. Locasale J.W. Metabolic landscape of the tumor microenvironment at single cell resolution Nat. Commun. 10 2019 3763 31434891
24 Jia Q. Wu W. Wang Y. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer Nat. Commun. 9 2018 5361 30560866
25 Mariathasan S. Turley S.J. Nickles D. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells Nature 554 2018 544 548 29443960
26 Gide T.N. Quek C. Menzies A.M. Wilmott, distinct immune cell populations define response to anti-PD-1 monotherapy and anti-PD-1/anti-CTLA-4 combined therapy Cancer Cell 35 2019 238 255.e6 30753825
27 Stuart T. Butler A. Hoffman P. Comprehensive integration of single-cell data Cell 177 2019 1888 1902.e21 31178118
28 Kröger H. Donner I. Skiello G. Influence of a new virostatic compound on the induction of enzymes in rat liver Arzneimittelforschung 25 1975 1426 1429 24
29 Morris T.J. Butcher L.M. Feber A. ChAMP: 450k chip analysis methylation pipeline Bioinformatics 30 2014 428 430 24336642
30 Shen S. Wang G. Shi Q. Seven-CpG-based prognostic signature coupled with gene expression predicts survival of oral squamous cell carcinoma Clin. Epigenetics 9 2017 88 28852427
31 Mayakonda A. Lin D. Assenov Y. Maftools: efficient and comprehensive analysis of somatic variants in cancer Genome Res. 28 2018 1747 1756 30341162
32 Hänzelmann S. Castelo R. Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data BMC Bioinf. 14 2013 7
33 Wilkerson M.D. Hayes D.N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking Bioinformatics 26 2010 1572 1573 20427518
34 Hu C. Li T. Xu Y. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data Nucleic Acids Res. 51 2023 D870 D876 36300619
35 Aibar S. González-Blas C.B. Moerman T. SCENIC: single-cell regulatory network inference and clustering Nat. Methods 14 2017 1083 1086 28991892
36 Liberzon A. Birger C. Thorvaldsdóttir H. The Molecular Signatures Database (MSigDB) hallmark gene set collection Cell Syst 1 2015 417 425 26771021
37 Jin S. Guerrero-Juarez C.F. Zhang L. Inference and analysis of cell-cell communication using CellChat Nat. Commun. 12 2021 1088 33597522
38 Love M.I. Huber W. Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 Genome Biol. 15 2014 550 25516281
39 Wu T. Hu E. Xu S. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data Innovation 2 2021 100141
40 Tommasini D. Fogel B.L. multiWGCNA: an R package for deep mining gene co-expression networks in multi-trait expression data BMC Bioinf. 24 2023 115
41 Geeleher P. Cox N. Huang R.S. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels PLoS One 9 2014 e107468
42 Yang W. Soares J. Greninger P. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells Nucleic Acids Res. 41 2013 D955 D961 23180760
43 Li J. Xing J. Yang Y. Adjuvant 131I-metuximab for hepatocellular carcinoma after liver resection: a randomised, controlled, multicentre, open-label, phase 2 trial Lancet Gastroenterol Hepatol 5 2020 548 560 32164877
44 Yau T. Park J.-W. Finn R.S. Nivolumab versus sorafenib in advanced hepatocellular carcinoma (CheckMate 459): a randomised, multicentre, open-label, phase 3 trial Lancet Oncol. 23 2022 77 90 34914889
45 El-Khoueiry A.B. Sangro B. Yau T. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial Lancet 389 2017 2492 2502 28434648
46 Yu M. Zeng W. Ouyang Y. ATP-exhausted nanocomplexes for intratumoral metabolic intervention and photoimmunotherapy Biomaterials 284 2022 121503
47 Hu Z. Qu G. Yu X. Acylglycerol kinase maintains metabolic state and immune responses of CD8+ T cells Cell Metab 30 2019 290 302.e5 31204281
48 Xie W. Chen B. Wen H. Biomimetic nanoplatform loading type I aggregation-induced emission photosensitizer and glutamine blockade to regulate nutrient partitioning for enhancing antitumor immunotherapy ACS Nano 16 2022 10742 10753 35830505
49 Sharabi A. Tsokos G.C. T cell metabolism: new insights in systemic lupus erythematosus pathogenesis and therapy Nat. Rev. Rheumatol. 16 2020 100 112 31949287
50 Ramapriyan R. Caetano M.S. Barsoumian H.B. Altered cancer metabolism in mechanisms of immunotherapy resistance Pharmacol. Ther. 195 2019 162 171 30439456
51 Sun R. Zhang Z. Bao R. Loss of SIRT5 promotes bile acid-induced immunosuppressive microenvironment and hepatocarcinogenesis J. Hepatol. 77 2022 453 466 35292350
52 Liu L. Liu Z. Gao J. CD8+ T cell trajectory subtypes decode tumor heterogeneity and provide treatment recommendations for hepatocellular carcinoma Front. Immunol. 13 2022 964190
53 Llovet J.M. De Baere T. Kulik L. Locoregional therapies in the era of molecular and immune treatments for hepatocellular carcinoma Nat. Rev. Gastroenterol. Hepatol. 18 2021 293 313 33510460
54 Singh P. Toom S. Avula A. The immune modulation effect of locoregional therapies and its potential synergy with immunotherapy in hepatocellular carcinoma J. Hepatocell. Carcinoma 7 2020 11 17 32104669
55 Galon J. Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies Nat. Rev. Drug Discov. 18 2019 197 218 30610226
56 Kalbasi A. Ribas A. Tumour-intrinsic resistance to immune checkpoint blockade Nat. Rev. Immunol. 20 2020 25 39 31570880
57 Pinyol R. Sia D. Llovet J.M. Immune exclusion-wnt/CTNNB1 class predicts resistance to immunotherapies in HCC Clin. Cancer Res. 25 2019 2021 2023 30617138
58 Long J. Wang A. Bai Y. Development and validation of a TP53-associated immune prognostic model for hepatocellular carcinoma EBioMedicine 42 2019 363 374 30885723
59 Kudo M. Finn R.S. Qin S. Lenvatinib versus sorafenib in first-line treatment of patients with unresectable hepatocellular carcinoma: a randomised phase 3 non-inferiority trial Lancet 391 2018 1163 1173 29433850
60 Bruix J. Qin S. Merle P. Regorafenib for patients with hepatocellular carcinoma who progressed on sorafenib treatment (RESORCE): a randomised, double-blind, placebo-controlled, phase 3 trial Lancet 389 2017 56 66 27932229
61 Granito A. Forgione A. Marinelli S. Experience with regorafenib in the treatment of hepatocellular carcinoma, Therap Adv. Gastroenterol. 14 2021 17562848211016959
