SPP1+ tumor-associated macrophages define a high-risk subgroup and inform personalized therapy in hepatocellular carcinoma

Introduction: Recently, contrary to attacking cancer cells, the tumor microenvironment (TME) with genomic stability and vulnerable nature has emerged as promising therapeutic targets in hepatocellular carcinoma (HCC). Within TME ecosystem, tumor-associated macrophages (TAMs) play a pivotal role in tumor evasion and progression of HCC. However, their clinical and therapeutic implications remain unexplored.

Methods: Utilizing a large-scale sc-RNA seq dataset, a landscape of HCC cellular ecosystem was depicted. Based on previous literature, an effectively differential TAMs subset classification was identified. Gene variations was extracted through trajectory analysis and then unsupervised clustering was conducted within RNA-seq data. Subsequently, survival analysis, specific pathway enrichment as well as hub regulatory network analysis were performed. Additionally, the immune cell infiltration and genomic variations were evaluated between clusters. Drug sensitivity and underlying therapeutic molecular were also explored. Through multiple immunofluorescence, our findings were verified.

Results: Herein, integrating single-cell RNA sequencing (scRNA-seq) and bulk RNA-seq data, we established a novel TAM classification system based on mutually exclusive SPP1 and FOLR2 signatures. According to the TAM trajectory genes, unsupervised clustering stratified HCC into three distinct clusters. Cluster 3 (C3), which is characterized by metabolic dysregulation and immunosuppressive TME, exhibited the poorest prognosis among the three groups. Hub network analysis of C3 further indicated its characteristic dysregulation of liver-specific metabolism. SPP1 was identified as a key signature of C3, which contributed to suppressing the infiltration of CD8⁺ T cells. Therapeutic evaluation revealed that C3 were sensitive to chemotherapy and tyrosine kinase inhibitors, while those C1 and C2 were more suitable for immunotherapy. Drug screening identified potential therapeutic compounds for each cluster.

Conclusion: This study redefines the heterogeneity of TAMs beyond the M1/M2 paradigm, linking the TAMs trajectory genes to HCC patient stratification. SPP1 blockade emerged as a strategy for counteracting immunosuppression, and cluster-specific therapies may optimize the management of HCC.

1 Introduction

Liver cancer has the sixth highest incidence and is ranked as the third leading cause of malignant mortality worldwide (1, 2), with hepatocellular carcinoma (HCC) accounting for 85%–90% (3, 4). According to the patients’ tumor stage, liver function, and performance status, a wide range of treatment options, including organ transplantation, surgical resection, radiation, and transarterial and systemic therapies, are alternatives for HCC (4, 5). With the approval of new agents as first- and second-line therapies, in particular the combination of atezolizumab with bevacizumab (T+A) as frontline standard care, the survival outcomes are superior than ever before (6, 7). Nevertheless, the prognosis remains dismal. Owing to the discouraging objective response rate, therapy resistance, and the high relapse probability, the overall 5-year survival rate is limited to 18% (8–10). Moreover, HCC displays high complexity and heterogeneity. Patients in the same clinical stage can respond differently to the same therapy (5, 11). Although unprecedented efforts based on multi-omics have been made for the molecular stratification of HCC, there is still a lack of a simple classification deriving from a single variable or molecular biomarker for effective clinical practice.

Currently, contrary to directly attacking cancer cells, the tumor microenvironment (TME) with genomic stability and vulnerable nature has increasingly emerged as a promising therapeutic target (12, 13). Among the cellular ecosystems of solid malignancy, tumor-associated macrophages (TAMs) play a pivotal role in governing cellular and molecular interactions, sustaining hallmarks of cancer, shaping TMEs, and eventually mapping onto clinical outcomes, as elucidated in studies (14, 15). Although the conventional M1/M2 macrophage classification simply and effectively defined the function and differentiation state of macrophages in vitro (13, 16), accumulating evidence based on single-cell RNA sequencing (scRNA-Seq) has challenged the applicability of this dichotomy for complex macrophages in vivo, particularly the TAMs (17, 18). Although multiple high resolutions of TAM subsets have been defined by scRNA-seq, the question of how to classify TAMs based on a relevant yet simple variable remains unanswered. It is well known that SPP1 encodes osteopontin (OPN), a phosphorylated glycoprotein regarded as a key component for tumor cell evolution and microenvironment reprogramming (19). Recently, SPP1+ macrophages have been recognized in several types of tumors and tend to present the malignant polarity of TAMs. The intratumoral cellular programs of SPP1+ TAMs, including promoting angiogenesis, enhancing tumor cell invasion, and resisting immune checkpoint blockade (ICB) therapy, have been reported (13, 14, 20–22). Nevertheless, how to define the clinical phenotypes and guide population-oriented therapy for HCC using SPP1+ TAM signatures has received little attention. There is an urgent need for the simple stratification of patients with HCC and the promotion of subset-targeted interventions.

scRNA-seq is a powerful technology for characterizing the heterogeneity of complex biological systems (23, 24). As a counterpart, the traditional RNA-seq approach, which considers the tumor but not the cell as a unit, obtains the average gene expression and individual phenotype (25). Trajectory analysis is widely used for cell differentiation inference and transcriptome dynamic process decoding (26). In particular, compared with superficial differential expression gene analysis that is conducted at a cluster resolution, trajectory analysis provides an effective way to explore intricate cell-to-cell variations at a single-cell resolution (13, 27). Herein, we depict a landscape of the cellular ecosystem of HCC utilizing a large-scale scRNA-seq dataset and identify a novel differential TAM subset classification for SPP1+ and FOLR2+ macrophages. Through trajectory analysis, we extract the underlying gene variations and define the TAM-based molecular classification of HCC, termed subgroups C1, C2, and C3, within the RNA-seq data. Furthermore, we explore the survival analysis, biological characteristics, hub regulatory gene network, and the genome variations between these three subgroups. We also evaluate the drug sensitivity and identify the underlying therapeutic molecules for each cluster. Accordingly, subtype-specific therapeutic strategies are proposed.

2 Materials and methods

2.1 Data acquisition and processing

A total of six independent cohorts were collected and processed in this study. The HCC scRNA-seq dataset GSE149614 was downloaded from the Gene Expression Omnibus (GEO). Excluding eight normal liver tissues (NLTs), we retained 10 primary tumor (PT), two portal vein tumor thrombus (PVTT), and one metastatic lymph node (MLN) samples (28). The Cancer Genome Atlas liver cancer cohort (TCGA-LIHC) and corresponding phenotype information, including 335 tumor samples and 89 normal samples, were retrieved from TCGA (https://xenabrowser.net/). A total of 221 RNA expression profiles of HCC from the GSE14520 cohort and 95 profiles from the GSE76427 cohort were obtained from the GEO database (29, 30). The International Cancer Genome Consortium—Liver Cancer–RIKEN Japan (ICGC-LIRI) and related phenotype information were downloaded from HCCDB v2.0 (http://lifeome.net:809/#/home) (31). In addition, the paraffin-embedded liver cancer tissue samples during general surgery were collected in the Mengchao Hepatobiliary Hospital of Fujian Medical University from January 2024 to December 2024. The study was conducted in accordance with the Declaration of Helsinki, and the study protocol was approved by the Institutional Review Board of the Mengchao Hepatobiliary Hospital of Fujian Medical University on November 3, 2023.

2.2 Single-cell RNA sequencing

The Seurat (version 4.4.0) package was used for processing and further analysis of data (32). Single-cell gene expression profiles were filtered using the criteria of a minimum threshold of 200 genes and a maximum threshold of 20% mitochondria genes per cell. The top 2,000 high variant genes were identified and scaled with the ScaleData function. To eliminate batch effects, all cells from the PT, PVTT, and MLN were integrated using Harmony (33, 34). The parameter setting in the Harmony algorithm was as follows: group.by.vars = “orig.ident”, max.iter.harmony = 10, lamba = 1. Using the RunPCA, FindNeighbor, and FindClusters functions, principal component analysis (PCA) linear dimensionality reduction and cluster visualization were performed. Cell clusters were annotated based on the cell lineage-specific genes (20, 28, 35), such as CD3D, CD3E, and CD3G for T/natural killer (NK) cells; PECAM1, DCN, and TM4SF1 for endothelial cells; ACTA2, DCN, and COL1A2 for fibroblasts; MS4A1, CD79A, and MZB1 for B/plasma cells; LYZ, CD14, and CD68 for myeloid cells; and ALB, SERPINA1, and KRT8 for hepatocytes. Sub-cluster analysis of myeloid cells was conducted with the SingleR package. According to previous documents, two different TAM clusters were identified based on specific cell type genes, including SPP1+ TAMs (SPP1, TREM2, FABP5, and NUPR1) and FOLR2+ TAMs (CD163, FOLR2, C1QB, and SEPP1) (13). These annotations were confirmed using the random forest algorithm. The “FindMarkers” function in Seurat was used to determine the differentially expressed genes (DEGs) between the two clusters. The non-parametric Wilcoxon rank-sum test was used to obtain p-values for comparison. The Monocle2 package was utilized to conduct pseudotime trajectory analysis of the TAMs for cell-to-cell variant demonstration (36). Genes along the trajectory were obtained and were enrolled in subsequent patient classification.

2.3 Identification of TAM-related subtypes

Non-negative matrix factorization (NMF) is an algorithm that can reduce high-dimensional datasets of tens of thousands of genes into a handful of metagenes, which are biologically easier to interpret (37). NMF has been widely used in various fields such as image analysis, speech recognition, auto signal processing, and bioinformatics (38). Using univariate Cox regression, prognosis-related genes were filtered from the trajectory genes. Subsequently, the consensus matrix of the filtered genes and clustering were constructed with the NMF algorithm. The parameter settings in the NMF algorithm were as follows: factorization ranks = 2–10, methods = “lee,” number of runs = 100 (38). The optimal rank was determined according to the cophenetic coefficient, before which a sharp decrease was observed (38). The silhouette graphic was used to qualify the clustering robustness. Specifically, individual silhouette statistics ranges from −1 to +1, where a high value indicates that the object is well matched to its own cluster and poorly matched to neighboring clusters (39).

2.4 Exploration of cluster-specific biological characteristics

Analysis of the DEGs between groups was performed using the DESeq2 package. The criteria log2 |foldchange| > 2 and p < 0.01 were used to filter genes. With the sorted expression profile, gene set enrichment analysis (GSEA) was conducted to elucidate the biological characteristics of distinct groups. The C2.cp.kegg.v5.2.symbols.gmt from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database was selected as the reference gene set. Hub gene-associated functional pathway enrichment was performed using the Gene Ontology (GO) pathways.

To evaluate the relative abundance of SPP1+ and FOLR2+ macrophages within the tumor and adjacent normal tissues, the DEGs [filtered using log2 (foldchange) > 0.585 and p < 0.01] in these two cell types were used as signature gene sets. Thereafter, gene set variation analysis (GSVA) (40) was applied to calculate the GSVA scores of the SPP1+ and FOLR2+ macrophages in TCGA and other independent datasets. Moreover, to explore the specific functional pathways of distinct macrophage clusters, the GSEA algorithm was performed using the R package clusterProfiler. The DEGs were determined using the “FindMarkers” function in the R package Seurat as described above.

2.5 Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) is commonly used to reveal the patterns of gene expression, summarizing the interconnections between modules and clinical traits, as well as identifying the candidate biomarkers or therapeutic targets (41, 42). In this study, cluster-associated regulated genes were investigated using the WGCNA package. Specifically, genes with the top 25% variance were selected and outlier samples were removed. By choosing an optimal soft threshold β (β = 5), the correlation matrix was converted into an adjacent matrix, and a topological overlap matrix (TOM) was subsequently formed. Utilizing average linkage hierarchical clustering, the 1 − TOM dissimilarity metric categorized genes with similar expression into gene modules. With the dynamic tree cut function, gene modules with diverse colors were determined. Finally, based on the relationship between modules and clinical traits, cluster-related modules with the tightest correlation were identified and signature genes were accordingly extracted.

2.6 Hub gene screening and visualization

The STRING website (https://cn.string-db.org/) is widely used for exploring protein interaction networks. Overlapping genes between the signature genes and the DEGs [filtered using log2 |foldchange| > 2 and p < 0.01] were submitted to the STRING website, and then the relative interaction network was retrieved, which was further imported into the Cytoscape software for visualization. The maximal clique centrality (MCC) algorithm in the cytohubba plug-in was used to obtain the top 10 ranked nodes as hub genes. Another plug-in, molecular complex detection (MCODE), was used to screen the hub modules (degree = 2, node score cutoff = 0.2, k-core = 2, max.depth = 100). A Venn diagram was used to capture an overlap of the hub nodes that were determined by the two algorithms mentioned above.

2.7 Immune cell infiltration analysis

xCell is a gene signature-based tool that utilizes a large compendium of publicly available transcriptomic data to infer the enrichment scores of 64 immune and stromal cell types (43). Using the xCell algorithm, the abundance of the adaptive immunity cells within each sample was quantified. Comparative analyses were performed to assess differential immune cell infiltration patterns across distinct clusters. Correlation analyses were conducted to investigate the relationship between the abundance of the adaptive immune cells and the expression levels of the cluster-specific signature genes.

2.8 Multiplexed immunofluorescence staining

Tumor tissues were first fixed in 10% formalin, embedded in paraffin, and then serially sectioned to 4-µm thickness. The following primary antibodies (all from ServiceBio, Wuhan, China) were used: CD8 (dilution 1:2,000), CD68 (dilution 1:5,000), and SPP1 (dilution 1:5,000). Subsequently, the samples were incubated with primary antibodies, followed by secondary antibodies. Double staining of CD8 with SPP1 and CD68 with SPP1 was mutually conducted. 4′,6-Diamidino-2-phenylindole (DAPI) was used for visualization of the cell nuclei (ServiceBio, Wuhan, China). A Nikon ECLIPSE C1 microscope was used for all imaging.

2.9 Somatic mutation analysis between clusters

The maftools package was utilized to identify and depict the top 10 somatic mutations among different clusters, including single nucleotide polymorphisms (SNPs), insertions and deletions (INDELs), the tumor mutation burden (TMB), and the mutation frequency (44).

2.10 Drug response prediction and potential therapeutic agent identification

Tumor Immune Dysfunction and Exclusion (TIDE) is commonly used to predict the immune-escape probability to comprehensively evaluate the T-cell infiltration and the T-cell function status. The TIDE score is available through the online algorithm (http://tide.dfci.harvard.edu/). In general, higher TIDE scores are associated with a poorer immune checkpoint inhibitor (ICI) therapeutic effect (45). The oncoPredict is an R package used to predict drug response and potential biomarkers based on cell line screening data (46). Importing the training matrices (GDSC and CTRP) from the website (https://osf.io/c6tfx), the sensitivity of the different clusters to chemotherapeutic or targeted drugs was predicted. The connectivity map (CMap) is a powerful tool based on pattern-matching algorithms for the identification of potential therapeutic compounds for specific populations (http://clue.io/) (47). The overlapping genes between the signature genes and the DEGs [filtered using log2 |foldchange| > 2 and p < 0.01) were selected and then the similarity of their expression profiles compared using the CMap database. Finally, the top 10 compounds from each cluster were obtained according to the ascending negative scores.

2.11 Statistical analysis

The Kruskal–Wallis test was used to compare the DEGs among the three groups. The survival probability between clusters was compared with a log-rank test. The Kaplan–Meier and Cox regression analyses were performed with the survival R package. Independent risk factors were identified using the multivariate Cox proportional hazards regression model. All statistical analyses were performed with R version 4.2.2.

3 Results

3.1 Single-cell transcriptomic atlas and myeloid cell landscape of HCC

To depict the landscape of the global microenvironment of HCC, single-cell transcriptomic analysis was performed on tumor sections from 10 patients, including 10 PTs, two PVTTs, and one MLN (Supplementary Table S1). After log-normalization and dimensionality reduction, a total of 43,228 cell transcriptomes were classified into 23 clusters (Figure 1A). According to the respective canonical markers, six main cell types—hepatocytes, T/NK cells, B/plasma cells, myeloid cells, fibroblasts, and endothelial cells—were defined (Figures 1B, C; Supplementary Figures S1A–D; Table S2). To further investigate the TAMs, a total of 8,038 myeloid cells were filtered out and categorized into three subpopulations: macrophages, monocytes, and dendritic cells (DCs) (Supplementary Figure S1F). Significantly, beyond the conventional M1/M2 macrophage dichotomy, two distinct macrophage phenotypes characterized by SPP1+ macrophage and FOLR2+ macrophage signatures were identified at the single-cell level. The t-distributed stochastic neighbor embedding (tSNE) map was used to delineate the segregated macrophage clusters, composed of the SPP1+ TAM and FOLR2+ TAM subsets, which exhibited mutually exclusive and DEG signatures (Figures 1D, E; Supplementary Figures S2F–H). These subset annotations were confirmed by analysis of the DEGs (Figure 1F) and the random forest algorithm (Supplementary Figure S2I). SPP1 emerged as one of the top upregulated genes in SPP1+ TAMs, but was significantly downregulated in FOLR2+ TAMs. Tumor samples were histologically accessed in order to gain insights into the distribution of SPP1+ macrophages in our cohort. Specifically, both SPP1+ microphages (with co-localization of CD68 and SPP1) and SPP1− macrophages (CD68⁺ only) were observed (Figure 1G).

Figure 1

Figure 1. Single-cell transcriptomic atlas of the global microenvironment and myeloid cell landscape of HCC. (A, B) UMAP plot of a total of 23 clusters (A) and six main cell types (B) from three tissue sites of 10 patients with HCC. (C) Dot plot showing the percentage of expressed cells and average expression levels of the canonical markers among six major cell types. Circle sizes represent the percentage of cells within a cluster expressing a gene. Color represents the average expression of each gene. (D, E) tSNE plots (D) and the mutually exclusive expression of the SPP1 and FOLR2 signatures (E) within myeloid cells. (F) Top 10 upregulated and downregulated genes of the two clusters. (G) mIHC staining determining the distribution of SPP1+ macrophages. (H, I) Signature scores of SPP1+ macrophages (H) and FOLR2+ macrophages (I) between tumor and adjacent normal tissues. HCC, hepatocellular carcinoma; UMAP, uniform manifold approximation and projection; tSNE, t-distributed stochastic neighbor embedding; mIHC, multiplex immunohistochemistry.

3.2 Trajectory analysis revealed TAM variations

Using the GSVA algorithm, SPP1+ macrophages were found to be significantly increased in tumor tissues, while FOLR2+ macrophages were enriched in adjacent normal tissues (Figures 1H, I; Supplementary Figures S2A, B). It was further observed that higher expression levels of SPP1 were remarkably associated with shorter overall survival (OS) in independent HCC cohorts (Figures 2A–D). Distinct functions of the two subtypes were also revealed. SPP1+ macrophages highly expressed metalloproteinases (e.g., MMP12), macrophage migration inhibitory factor (e.g., MIF), and chemokines (e.g., CXCL2 and CXCL3). On the other hand, FOLR2+ macrophages showed higher expression of CD5L and the metabolism-associated genes (e.g., APOA2, APOA3, and AHSG) (Figure 1F). Functional analysis showed that the chemokine signaling pathway, cytokine–cytokine receptor interaction, glycolysis/gluconeogenesis, the IL-17 signaling pathway, and the PPAR signaling pathway were enriched in SPP1+ macrophages, while complement and coagulation cascades, cholesterol metabolism, and metabolism of xenobiotics by cytochrome P450 were highly enriched in FOLR2+ macrophages (Supplementary Figure S2C).

Figure 2

Figure 2. Trajectory analysis revealing variations in TAMs. (A–C) Plots of the Kaplan–Meier overall survival curves from TCGA-LIHC (A), GSE14520 (B), and ICGC-LIRI (C) cohorts, grouped with high and low expression of SPP1. (D) Univariate analyses of SPP1 within the different datasets. (E) Pseudotime trajectory showing the distribution of macrophages. (F) The SPP1 signature accumulated with the progress of pseudotime. (G) FOLR2 signature enriched in the early trajectory stage. (H) Heatmap displaying the genes along with a dynamic trend. (I) Peak map showing the distribution of the SPP1+ and FOLR2+ TAMs during the pseudotime transition. (J) Identification the SPP1+ TAM-related population. (J) Cophenetic curve revealing 3 as the optimal rank. (K) Consensus matrix map of the NMF clustering of the HCC cases from the TCGA-LIHC cohort. (L) Silhouette plot validating the robustness of the NMF cluster. TAMs, tumor-associated macrophages; NMF, non-negative matrix factorization; HCC, hepatocellular carcinoma; TCGA-LIHC, The Cancer Genome Atlas—Liver Hepatocellular Carcinoma; ICGC-LIRI, International Cancer Genome Consortium—Liver Cancer–RIKEN Japan.

For in-depth exploration of the phenotypic dynamics of the macrophages and their correlated clinical relevance, the Monocle algorithm was applied, considering cells, not clusters, as a unit of statistical replication (27). Notably, it was found that the SPP1 expression progressively accumulated along a pseudotime trajectory, while the FOLR2 signature was predominantly enriched in the early stage (Figures 2E–I). Within the trajectory, molecules with similar dynamic trends were gathered and displayed in a pseudotime heatmap. The genes along the trajectory were extracted and two coherently regulated gene modules were identified (Figure 2H). Notably, SPP1 was determined as one of the main dynamic genes (Figure 2H; Supplementary Table S3; the heatmap shows the top 50 trajectory genes).

3.3 Identification of the SPP1+ TAM-related population

For more in-depth insights into the clinical phenotype changes along with dynamic shifts of the macrophage subtypes, a non-negative matrix was constructed according to the prognosis-associated genes filtered from the trajectory genes (Supplementary Table S4). For the sharpest decrease in coefficient observed before 4 in the cophenetic arrangement, 3 was determined as the optimal rank (Figure 2J). Unsupervised clustering analysis of a total of 365 tumor patients within TCGA-LIHC determined three groups: 232 cases in cluster 1 (C1), 45 cases in cluster 2 (C2), and 88 cases in cluster 3 (C3). Significant heterogeneity of the transcriptome was observed between clusters (Figure 2K). The silhouette graphic further provided proof the robustness of the stratification (Figure 2L). In addition, survival analysis displayed that C1 yielded favorable OS, while C3 tended to have the poorest prognosis. Using a multivariate Cox regression model, a C3 classification was determined to be an independent risk factor for shorter survival, whereas C1 was deemed as a protective factor for better outcomes (Figures 3A, B). Similarly, the bar graph illustrated that C1 and C2 predominantly comprised patients in the early stage (I/II), while the C3 stratification was significantly associated with patients in the advanced stage (III/IV) (p < 0.05; Figure 3C), suggesting potential tumor progression from the C1 to the C3 phenotype. Moreover, the underlying biological characteristics were delineated between clusters using GSEA. Consistent with the FOLR2+ macrophage signature, C1 exhibited activation of the metabolic pathways including fatty acid degradation (FAD), amino acid metabolism (glycine, serine, and threonine metabolism), and metabolism of xenobiotics by cytochrome P450 (Figure 3D). Conversely, these pathways were significantly suppressed in C3. Notably, previous research revealed FAD to be inhibited in multiple cancers and correlated with increased fatty contents in tumors, which impaired the metabolism of fatty acids (FAs) and promoted the pro-tumoral phenotype polarization of TAMs. Although the amino acid pathway is typically aberrantly activated during oncogenesis, its opposite dysregulation may suppress T-cell proliferation and antitumor immunity. For C3, several significantly upregulated pathways, including ECM–receptor interaction and the IL-17 signaling pathway, were noticeably enriched, similarly to the functional profile of SPP1+ macrophages (Supplementary Figure S2C). Overall, C1 was characterized as a metabolism-enhanced, FOLR2+ macrophage-related HCC subtype, while C3 was characterized as a metabolism-dysregulated, SPP1+ macrophage-associated HCC subtype. Metabolic alterations and TME remodeling underlie the observed clinical heterogeneity.

Figure 3

Figure 3. Exploration of the SPP1+ TAM-related population. (A) Kaplan–Meier curves of the overall survival of three subgroups in TCGA-LIHC. (B) Forest plot of the survival-related factors. (C) Bar plot displaying the correlation between HCC stage and cluster stratification. (D) GSEA showing the enriched signaling pathways of each cluster. (E) The Dendrogram displaying the gene modules with diverse colors. (F) Correlation analysis of module eigengenes and molecular phenotypes. TCGA-LIHC, The Cancer Genome Atlas—Liver Hepatocellular Carcinoma; HCC, hepatocellular carcinoma; GSEA, gene set enrichment analysis.

3.4 Investigation of the hub molecule network embedded in the SPP1+ TAM-related cluster

To explore the molecular characteristics of the SPP1+ macrophage-associated subtype, WGCNA was performed. After filtering qualified genes, stratifying all samples, defining an optimal soft threshold β (β = 5), and transforming the matrix, a total of seven co-expression modules were identified (Figure 3E; Supplementary Figures S3A, B). The correlations between the gene modules and the clinical traits were determined using a module–trait relationship heatmap (Figure 3F; Supplementary Figure S3C). The turquoise module exhibited the strongest association with C3. The correlation coefficient between module membership and gene significance validated the robustness of the characteristic genes (Figure 4A). Furthermore, the hub genes within the turquoise module were used to construct a gene interaction network using the STRING website (Figure 4B; Supplementary Table S5). Using the Maximal Clique Centrality (MCC) and Molecular Complex Detection (MCODE) algorithms, hub genes and the core module of the interaction network were identified. Remarkably, a significant overlap was observed between the nodes of the core module and the top 10 hub genes (Figures 4C, D; Supplementary Table S6), indicating their critical role in this network.

Figure 4

Figure 4. Investigation of the hub molecular network of the SPP1+ TAM-related cluster. (A) Scatter plot showing the module membership and gene significance within the turquoise (C3-related) module. (B) Interaction network analysis of the genes in the turquoise module using Cytoscape software. (C) Venn diagram displaying the overlap of the hub genes and hub module nodes. (D) Interaction network of the five overlapping hub genes. (E) Correlation of the hub gene expression with overall survival (OS) in The Cancer Genome Atlas—Liver Hepatocellular Carcinoma (TCGA-LIHC). The blue line designates samples with lowly expressed genes, while the yellow line indicates samples with highly expressed genes. (F) Bar plot showing the Gene Ontology (GO) enrichment pathways of the hub genes in the C3-related module. (G) The most important hub genes and their corresponding GO pathways. TAM, tumor-associated macrophage; MCODE, molecular complex detection.

Furthermore, the expression patterns of these hub genes were determined in tumor versus paraneoplastic tissues, as well as across C3 and the other subgroups. Specifically, UGT1A10 exhibited a markedly elevated expression in tumor tissues, while the remaining nine hub genes (ADH1A, ADH1B, ADH1C, ADH4, SLC10A1, CYP2C9, GPT, HRG, and PON1) showed significantly reduced expression levels in cancerous tissues compared with their adjacent normal counterparts (Supplementary Figure S3E). A similar trend in expression was also observed when comparing the C3 subgroup with the other subgroups, hinting at the differential expression patterns of these genes in tumor context (Supplementary Figure S4A). Based on the expression patterns in cancerous and non-cancerous tissues, the survival analysis revealed that the decreased expression of nine hub genes—ADH1A, ADH1B, ADH1C, ADH4, SLC10A1, CYP2C9, GPT, HRG, and PON1—was significantly associated with worse prognostic outcomes. The elevated expression of UGT1A10 was correlated with shorter survival time (Figure 4E). Data from the Human Protein Atlas revealed that most of the hub genes (i.e., ADH1A, ADH1B, ADH1C, ADH4, CYP2C9, GPT, HRG, and PON1) belong to the liver-specific enzyme family that is involved in metabolizing various xenobiotic compounds, including alcohol, retinol aliphatic alcohols, hydroxysteroids, and lipid peroxidation products. Pathway analysis also confirmed the involvement of the hub genes in multiple metabolic alterations (Figures 4F, G). The downregulation of these hub genes may underlie the characteristic metabolic dysregulation of C3.

Moreover, for C3, which showed a marked enrichment of immune- and inflammation-associated pathways (Figure 3D), the potential role of the hub genes in the regulation of immune response was analyzed. Notably, decreased expression of the hub genes was significantly associated with the decreased infiltration of adaptive immune cells, particularly CD8⁺ naive T cells and activated NK cells (Supplementary Figure S4B).

3.5 Determination of the cluster-driven signature of HCC population stratification

To investigate the molecular signatures that drive the HCC population stratification, a heatmap was constructed based on the cluster-related module genes (Figure 5A). Exhibiting diverse expression patterns across the HCC groups, SPP1 emerged as a predominant cluster-driven signature. A remarkably elevated expression of SPP1 was observed in C3 (Figure 5B). Based on the downregulated hub genes within C3, a negative cluster score (−C3 signature score) was calculated to characterize this subgroup. Notably, the expression level of SPP1 showed a strong correlation with the −C3 signature score (Figure 5C), further confirming SPP1 as a key signature of the C3 subgroup. To investigate the mechanisms underlying the association between the C3 population and poor prognosis, the infiltration of adaptive immune cells was assessed using the xCell algorithm. Intriguingly, the C3 subgroup exhibited a significant reduction in CD8⁺ naive T-cell infiltration compared with the other subgroups (Figures 5D–F). Furthermore, utilizing the TIMER 2.0 database, it was demonstrated that the expression of SPP1 was inversely correlated with the infiltration of both CD8⁺ naive T cells and total CD8⁺ central memory T cells in the HCC environment (Supplementary Figures S4C, D). Consistently, multiplex immunohistochemistry (mIHC) staining showed that SPP1 tended to localize at the tumor boundary, whereas CD8⁺ T cells were likely to localize outside of the tumor, with less infiltration in tumor tissues (Figure 5G), hinting at a potential role of SPP1 in shaping the immunosuppressive microenvironment of HCC and at a promising strategy to counteract immunosuppression.

Figure 5

Figure 5. Identification of the cluster-driven signature in C3 hepatocellular carcinoma (HCC) stratification. (A) Heatmap depicting the driven signatures of the HCC clusters, with SPP1 emerging as a predominant cluster-driven signature. (B) Violin plot of the SPP1 expression levels between clusters. (C) Correlation between the expression of SPP1 and the negative cluster score (−C3 signature score). (D–F) Bar plots showing the infiltration levels of CD8⁺ naive T cells (D), CD8⁺ effector memory T cells (E), and CD8⁺ central memory T cells (F) across three clusters. (G) Multiplex immunohistochemistry (mIHC) staining revealing the characteristic localization of SPP1 and CD8 in HCC, with SPP1 tending to suppress the infiltration of CD8⁺ T cells.

3.6 Identification of subtype-specific therapeutic strategies and promising treatment targets

As distinct molecular phenotypes enable personalized treatment and tailored clinical management for specific HCC subgroups, potential precision strategies for macrophage-associated subtypes were further explored. As is known, transarterial chemoembolization (TACE) remains the preferred treatment for unresectable intermediate-to-advanced HCC. Using the oncoPredict algorithm, the sensitivity of the different clusters for conventional chemotherapeutic drugs was evaluated. As shown in Figure 6A, with the relatively lower half-maximal inhibitory concentration (IC₅₀), C3 was more sensitive to most chemotherapeutic agents, including 5-fluorouracil, cytarabine, camptothecin, docetaxel, gemcitabine, and epirubicin. The results are in line with the aforesaid stage stratification (Figure 3C), wherein C3 was dominated by a population in the advanced tumor stage (III/IV). The sensitivity of the different groups to targeted therapy agents was also determined. Significantly, patients in C3 exhibited superior response to multiple kinase inhibitors, including sorafenib, alpelisib, and gefitinib, compared with those in the other groups. Given the critical role of somatic mutation in shaping cancer phenotypes and the therapeutic response, the mutational landscapes were characterized using the maftools algorithm. Consistent with the established literature, TP53, CTNNB1, and MUC16 were identified as the most frequently mutated genes across all subgroups. C3 demonstrated a distinctive mutational profile, with predominant TP53 mutations (21%, 27%, and 47% positivity, respectively), in contrast to the CTNNB1 predominance observed in clusters 1 and 2 (29%, 34%, and 14% positivity, respectively) (Figure 6B), hinting at the instability of the genome and the dysregulation of the cell cycle in C3. In addition, to further delineate the populations suitable for immunotherapy, the TIDE scores were calculated. A significant increase in the TIDE score was observed in C3, suggesting a poorer response to ICIs. The other two subgroups maintained comparatively lower TIDE scores (Figure 6C), implicating a potential group for ICIs. Moreover, using the CMap, potential therapeutic compounds were identified for each cluster (Figure 6D; Supplementary Figures S4A, C), and correspondingly, the underlying regulatory mechanisms were also elucidated (Supplementary Figures S4B, D). Collectively, according to the distinct cluster stratification, personalized treatment strategies would be available and ultimately improve clinical practice.

Figure 6

Figure 6. Evaluation of clinical treatment and identification of promising therapeutic targets. (A) Box plot of the estimated IC₅₀ values of chemotherapy and targeted therapy drugs. (B) Waterfall plots showing the distribution of the somatic mutation genes with the top 10 highest mutation frequencies. (C) TIDE algorithm for the prediction of the immunotherapy response between clusters. (D) Potential therapeutic compounds for the C3 group. IC₅₀, half-maximal inhibitory concentration; TIDE, Tumor Immune Dysfunction and Exclusion.

4 Discussion

Although multiple molecular classifications have improved our understanding of the heterogeneity of HCC, the complicated classifier and the reproducible difficulty still pose a substantial challenge in the translation of the subtypes into clinical practice (11, 48, 49). In contrast to previous classifications that focused on comprehensive genomic or transcriptomic analysis, herein, we were concerned on the dynamic evolution of TAMs, potential heterogeneous regulators (14, 15, 50). We specifically identified a novel stratification of TAMs, proposed a simple classification of HCC, and developed personalized therapy correspondingly.

Currently, TMEs with genomic stability and vulnerable nature have increasingly emerged as promising therapeutic targets (12, 13). Multiple associated attempts, including T cells with cytolytic functions or combination therapies that inhibit compensatory signaling pathways, have been developed (9, 51, 52). In this study, using scRNA-seq analysis, we revealed the main cellular components of the TMEs within HCC. We typically focused on the TAMs, pivotal components revealed to govern the cellular and molecular interactions and to sustain hallmarks of cancer, ultimately mapping onto clinical outcomes (14, 15, 50). Similarly to recent single-cell transcription analysis that abandoned the traditional pro-inflammatory M1 and anti-inflammatory M2 dichotomy (50, 53), we failed to clearly distinguish the M1 and M2 macrophages using known markers, such as CD86 (M1) and CD163 (M2). Notably, based on previous knowledge (13, 14), we observed a roughly exclusive expression of the SPP1 and FOLR2 genes within the macrophage clusters. With regard to the complexity of the cellular ecosystem and the instability of the sub-clusters derived from the dimensionality reduction, we innovatively annotated macrophages into SPP1+ and FOLR2+ TAMs according to simple but distinct variants. We validated the distribution of SPP1+ macrophages in our cohort through immunostaining. Diverse cell states were further defined by trajectory analysis, in which the frequency of FOLR2 became relatively diluted among TAMs by the expression of SPP1. We also observed significant enrichment of the SPP1+ macrophages in tumor tissues compared with adjacent normal tissues, where FOLR2+ macrophages predominated. As FOLR2+ macrophages were revealed to be embryonic-origin tissue-resident macrophages (TRMs, also known as Kupffer cell-like phenotypes in the liver) (50), it was implied that, during tumor progression, SPP1+ macrophages increased from the infiltration of circulating monocytes, accompanied with reduced FOLR2+ TRMs. For functional analysis, higher enrichment of the PPAR pathway was observed in SPP1+ TAMs, and increased expression of the downstream protein metalloproteinases (e.g., MMP12) of this pathway was also identified in SPP1+ macrophages, hinting at a potential function of the activation of the PPAR pathway in the terminal differentiation of SPP1+ TAMs. The dynamic evolution of TAMs may initiate intratumoral cellular program reshaping and heterogenetic clinical phenotype formation. The genes along the trajectory of TAMs provided revealing features that covaried between individuals, which enables a population-oriented molecular classification (54, 55).

Therefore, using the NMF algorithm, heterogeneous molecular clusters of HCC were obtained. Notably, the C3 group was associated with unfavorable prognosis and was regarded as an independent risk factor for worse outcomes. Several upregulated pathways, including the ECM–receptor interaction and the IL-17 signaling pathway, were simultaneously enriched in C3 and SPP1+ macrophages; therefore, C3 was assumed as an SPP1+ macrophage-associated HCC subtype. Significant metabolic alterations, including the downregulated FAD and amino acid metabolism (glycine, serine, and threonine metabolism), were observed in C3. These may be attributed to the reciprocal interactions between the dynamic TAM evolution and the metabolic reprogramming within the TMEs. Similar to de novo FA synthesis (56), the inhibition of FAD, as another energy resource, caused FA accumulation in the microenvironment, subsequently inducing dysregulated FA metabolism and promoting pro-tumoral TAM phenotype polarization (57, 58). Furthermore, it is well known that T cells in the TME catabolize lipids through mitochondrial fatty acid oxidation (FAO) to meet the energy demands under nutrient stress (59); however, current evidence demonstrates that unmodified cellular therapy products fail to sustain the bioenergetics in tumors (60). Within increasing FA concentrations due to the inhibition of FAD in the TME, we supposed that a compromised metabolic state of T cells with impaired antitumor effector function is formed (61, 62). Taken together, these mechanisms mobilized the immunosuppressive TME, facilitating tumor progression and proliferation. Although the amino acid pathways are typically upregulated during tumorigenesis, an inverse pattern of dysregulation was observed in the C3 group. This aberrant metabolic reprogramming resulted in altered energy homeostasis, which has been shown to suppress T-cell proliferation and attenuate the antitumor immune response (63, 64).

In addition, most of the hub genes in the C3-associated module belong to liver-specific enzymes and are predominantly enriched in metabolic alteration, which further confirmed the characteristic metabolic dysregulation in C3. Among these hub genes, typically, ADH1A is an enzyme that is involved in metabolizing various xenobiotic substrates (65, 66). Evidence suggests that the downregulation of ADH1A may facilitate transition from liver damage to hepatocarcinogenesis and exacerbate HCC progression upon exposure to xenobiotic compounds (67). Consistently, a reduction in the expression of ADH1A has been observed in hepatocarcinoma tissue compared with para-tumor tissues, and this pattern is also evident in C3 compared with C1 and C2. Similar influences have also been observed in other hub genes within C3, which conclusively contribute to the C3 phenotype formation.

Based on the downregulated hub genes, a negative cluster score (−C3 signature score) was defined for the C3 subgroup of HCC. For the characteristic gene driving the specific molecular classification, SPP1, derived from the overlap of the trajectory genes and the main module genes and significantly correlated with the −C3 signature score, was identified as the primary signature of C3. Also known as OPN, SPP1 is an integrin-binding glycoprotein that has been reported to be overexpressed in various tumors (68). The literature has elucidated the crucial role of SPP1 in remodeling the TME, including promoting migration and colony formation and facilitating M2-like polarization and immune cell suppression (69, 70). Our findings demonstrate the remarkable increase in the expression of SPP1 in HCC tissues, which was strongly correlated with a reduced CD8⁺ T-cell infiltration. The mIHC staining revealed that SPP1 tended to localize at the tumor boundary, whereas CD8⁺ T cells were likely to localize outside of the tumor, further confirming the immunosuppressive role of SPP1. Typically, among the three classification phenotypes, the expression of SPP1 was markedly elevated in the C3 compared with the C1 and C2 subtypes, which may partially explain the immunosuppressive microenvironment within C3. As ICI-based therapy has been established as the first-line treatment and its efficacy is highly dependent on TME (4), targeting SPP1 may be a promising strategy to regulate the TME and improve immunotherapy response.

As distinct molecular phenotypes uncovered features that covary between individuals, personalized treatment recommendations and tailored management are available for specific groups of HCC. Through the oncoPredict algorithm, the therapeutic efficacy of various treatments across the three subtypes was evaluated. Specifically, C3 exhibited greater sensitivity to most chemotherapeutic agents, including 5-fluorouracil, cytarabine, camptothecin, docetaxel, gemcitabine, and epirubicin. A plausible explanation for this observation may lie in the impaired metabolic pathways within C3, which imposed significant chemotherapeutic stress and influenced the subsequent treatment response (71). Sorafenib is the standard systemic therapy for advanced HCC (4). Notably, the superior efficacy of sorafenib was observed in C3 compared with the other subgroups, along with other kinase inhibitors including alpelisib and gefitinib. Previous studies have elucidated that sorafenib-mediated lipotoxicity contributes to its therapeutic efficacy, in which sorafenib acts as a direct LXR signaling activator and promotes lipogenesis and a toxic accumulation of FAs (11, 52). For patients in C3, the downregulated FAD pathway may have already predisposed them to lipotoxicity stress prior to treatment. The additional administration of sorafenib exacerbated this lipotoxic effect, which ultimately enhanced the therapeutic response in this group. Furthermore, consistent with data indicating that HCC with TP53 mutation is markedly associated with immunosuppression (67), it was confirmed that the C3 subgroup, which harbored the highest frequency of TP53 mutations, presented poorer response to immunotherapy. Given the lack of precision medicine in HCC, we identified novel potential therapeutic agents and their corresponding underlying mechanisms based on similar transcription patterns embedded in each cluster. Individualized therapy may be available, accordingly.

Several limitations of the present study need to be acknowledged. Primarily, all of the datasets enrolled in this study were retrospective. Secondly, we provide preliminary evidence that the SPP1+ and FOLR2+ TAMs, rather than the traditional M1/M2 dichotomy, represent a novel macrophage stratification at single-cell resolution. Adequately powered and well-designed studies are required to confirm these findings. In addition, although we revealed the tumor heterogeneity and stratified HCC into three molecular clusters, along with providing potential therapeutic recommendations for specific subgroups, the causality and the precise underlying mechanisms remain to be further elucidated. Population-oriented scRNA-seq analysis is expected for a better understanding and interpretation of our findings.

5 Conclusion

We established a novel TAM stratification system and classified HCC populations based on the TAM trajectory genes. SPP1 was identified as a key signature of the malignant cluster, which negatively impacted immune cell infiltration. Targeting SPP1 may improve the therapeutic efficacy, and personalized treatment strategies can be tailored based on specific patient stratification.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: The Cancer Genome Atlas portal (https://portal.gdc.cancer.gov/, accession number TCGA-LIHC), and Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/, accession number GSE149614, GSE14520 and GSE76427), and HCCDB v2.0 (http://lifeome.net:809/#/home, accession number ICGC-LIRI).

Ethics statement

The studies involving humans were approved by Institutional Review Board of the Mengchao Hepatobiliary Hospital of Fujian Medical University. The studies were conducted in accordance with the local legislation and institutional requirements. The human samples used in this study were acquired from primarily isolated as part of your previous study for which ethical approval was obtained. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and institutional requirements.

Author contributions

W-XX: Writing – original draft, Software, Writing – review & editing, Investigation, Formal analysis, Validation, Visualization, Data curation, Conceptualization. Y-MY: Visualization, Data curation, Writing – review & editing. J-LC: Validation, Methodology, Writing – review & editing, Data curation. X-YG: Methodology, Validation, Data curation, Writing – review & editing. CL: Methodology, Writing – review & editing. JL: Writing – review & editing, Methodology. L-BL: Visualization, Investigation, Funding acquisition, Formal analysis, Project administration, Resources, Software, Writing – review & editing. XC: Supervision, Conceptualization, Writing – review & editing, Funding acquisition, Writing – original draft, Resources, Project administration.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study were supported by Central Guidance on Local Science and Technology Development Fund of Fujian Province, No. 2023L3028; Natural Science Foundation of Fujian Province, No. 2021J01281; Science and Technology Innovation Joint Foundation of Fujian Province, No. 2024Y9421; Research Initiation Foundation of Mengchao Hepatobiliary Hospital of Fujian Medical University, No. QDZJ-2022-001; Health and Wellness Youth Scientific Research Project of Fujian Province, No. 2023QNA081.

Acknowledgments

The authors thank the staff of the Department of Oncology of Mengchao Hepatobiliary Hospital of Fujian Medical University, for their guidance and support.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1606195/full#supplementary-material

Abbreviations

HCC, hepatocellular carcinoma; TME, tumor microenvironment; ICB, immune checkpoint blockade; scRNA-seq, single-cell RNA sequencing; NLT, normal liver tissue; PT, primary tumor; PVTT, portal vein tumor thrombus; GEO, Gene Expression Omnibus; PCA, principal component analysis; NMF, non-negative matrix factorization; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; WGCNA, weighted gene co-expression network analysis; TOM, topological overlap matrix; MCC, maximal clique centrality; MCODE, molecular complex detection; DAPI, 4′,6-diamidino-2-phenylindole; SNP, single nucleotide polymorphism; INDEL, insertion and deletion; TMB, tumor mutation burden; TIDE, tumor immune dysfunction and exclusion; UMAP, uniform manifold approximation and projection; tSNE, t-distributed stochastic neighbor embedding; TAMs, tumor-associated macrophages; TCGA-LIHC, The Cancer Genome Atlas—Liver Hepatocellular Carcinoma; GSEA, gene set enrichment analysis; FDA, fatty acid degradation; ECM, extracellular matrix; mIHC, multiplex immunohistochemistry; TACE, transarterial chemoembolization; IC₅₀, half-maximal inhibitory concentration; ICIs, immune checkpoint inhibitors; CMap, connectivity map; TRMs, tissue-resident macrophages; FAs, fatty acids.

References

1. Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2024) 74:229–63. doi: 10.3322/caac.21834

PubMed Abstract | Crossref Full Text | Google Scholar

2. Siegel RL, Giaquinto AN, and Jemal A. Cancer statistics, 2024. CA Cancer J Clin. (2024) 74:12–49. doi: 10.3322/caac.21820

PubMed Abstract | Crossref Full Text | Google Scholar

3. Sartorius K, Sartorius B, Aldous C, Govender PS, and Madiba TE. Global and country underestimation of hepatocellular carcinoma (HCC) in 2012 and its implications. Cancer Epidemiol. (2015) 39:284–90. doi: 10.1016/j.canep.2015.04.006

PubMed Abstract | Crossref Full Text | Google Scholar

4. Vogel A, Meyer T, Sapisochin G, Salem R, and Saborowski A. Hepatocellular carcinoma. Lancet. (2022) 400:1345–62. doi: 10.1016/S0140-6736(22)01200-4

PubMed Abstract | Crossref Full Text | Google Scholar

5. Reig M, Forner A, Rimola J, Ferrer-Fabrega J, Burrel M, Garcia-Criado A, et al. BCLC strategy for prognosis prediction and treatment recommendation: The 2022 update. J Hepatol. (2022) 76:681–93. doi: 10.1016/j.jhep.2021.11.018

PubMed Abstract | Crossref Full Text | Google Scholar

6. Vogel A, Martinelli E, and Committee EG. Updated treatment recommendations for hepatocellular carcinoma (HCC) from the ESMO Clinical Practice Guidelines. Ann Oncol. (2021) 32:801–5. doi: 10.1016/j.annonc.2021.02.014

PubMed Abstract | Crossref Full Text | Google Scholar

7. Bruix J, Chan SL, Galle PR, Rimassa L, and Sangro B. Systemic treatment of hepatocellular carcinoma: An EASL position paper. J Hepatol. (2021) 75:960–74. doi: 10.1016/j.jhep.2021.07.004

PubMed Abstract | Crossref Full Text | Google Scholar

8. Ray K. Liver cancer: Nivolumab: checkmate for hepatocellular carcinoma? Nat Rev Gastroenterol Hepatol. (2017) 14:326. doi: 10.1038/nrgastro.2017.62

PubMed Abstract | Crossref Full Text | Google Scholar

9. Finn RS, Qin S, Ikeda M, Galle PR, Ducreux M, Kim TY, et al. Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. N Engl J Med. (2020) 382:1894–905. doi: 10.1056/NEJMoa1915745

PubMed Abstract | Crossref Full Text | Google Scholar

10. Craig AJ, von Felden J, Garcia-Lezana T, Sarcognato S, and Villanueva A. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. (2020) 17:139–52. doi: 10.1038/s41575-019-0229-4

PubMed Abstract | Crossref Full Text | Google Scholar

11. Li B, Li Y, Zhou H, Xu Y, Cao Y, Cheng C, et al. Multiomics identifies metabolic subtypes based on fatty acid degradation allocating personalized treatment in hepatocellular carcinoma. Hepatology. (2024) 79:289–306. doi: 10.1097/HEP.0000000000000553

PubMed Abstract | Crossref Full Text | Google Scholar

12. Xiao Y and Yu D. Tumor microenvironment as a therapeutic target in cancer. Pharmacol Ther. (2021) 221:107753. doi: 10.1016/j.pharmthera.2020.107753

PubMed Abstract | Crossref Full Text | Google Scholar

13. Bill R, Wirapati P, Messemaker M, Roh W, Zitti B, Duval F, et al. CXCL9:SPP1 macrophage polarity identifies a network of cellular programs that control human cancers. Science. (2023) 381:515–24. doi: 10.1126/science.ade2292

PubMed Abstract | Crossref Full Text | Google Scholar

14. Nalio Ramos R, Missolo-Koussou Y, Gerber-Ferder Y, Bromley CP, Bugatti M, Nunez NG, et al. Tissue-resident FOLR2(+) macrophages associate with CD8(+) T cell infiltration in human breast cancer. Cell. (2022) 185:1189–207 e25. doi: 10.1016/j.cell.2022.02.021

PubMed Abstract | Crossref Full Text | Google Scholar

15. Tan J, Fan W, Liu T, Zhu B, Liu Y, Wang S, et al. TREM2(+) macrophages suppress CD8(+) T-cell infiltration after transarterial chemoembolisation in hepatocellular carcinoma. J Hepatol. (2023) 79:126–40. doi: 10.1016/j.jhep.2023.02.032

PubMed Abstract | Crossref Full Text | Google Scholar

16. Pittet MJ, Michielin O, and Migliorini D. Clinical relevance of tumour-associated macrophages. Nat Rev Clin Oncol. (2022) 19:402–21. doi: 10.1038/s41571-022-00620-6

PubMed Abstract | Crossref Full Text | Google Scholar

17. Qian BZ and Pollard JW. Macrophage diversity enhances tumor progression and metastasis. Cell. (2010) 141:39–51. doi: 10.1016/j.cell.2010.03.014

PubMed Abstract | Crossref Full Text | Google Scholar

18. Locati M, Curtale G, and Mantovani A. Diversity, mechanisms, and significance of macrophage plasticity. Annu Rev Pathol. (2020) 15:123–47. doi: 10.1146/annurev-pathmechdis-012418-012718

PubMed Abstract | Crossref Full Text | Google Scholar

19. Rittling SR and Singh R. Osteopontin in immune-mediated diseases. J Dent Res. (2015) 94:1638–45. doi: 10.1177/0022034515605270

PubMed Abstract | Crossref Full Text | Google Scholar

20. Liu Y, Xun Z, Ma K, Liang S, Li X, Zhou S, et al. Identification of a tumour immune barrier in the HCC microenvironment that determines the efficacy of immunotherapy. J Hepatol. (2023) 78:770–82. doi: 10.1016/j.jhep.2023.01.011

PubMed Abstract | Crossref Full Text | Google Scholar

21. Qi J, Sun H, Zhang Y, Wang Z, Xun Z, Li Z, et al. Single-cell and spatial analysis reveal interaction of FAP(+) fibroblasts and SPP1(+) macrophages in colorectal cancer. Nat Commun. (2022) 13:1742. doi: 10.1038/s41467-022-29366-6

PubMed Abstract | Crossref Full Text | Google Scholar

22. Lyu A, Fan Z, Clark M, Lea A, Luong D, Setayesh A, et al. Evolution of myeloid-mediated immunotherapy resistance in prostate cancer. Nature. (2024) 637:1207–17. doi: 10.1038/s41586-024-08290-3

PubMed Abstract | Crossref Full Text | Google Scholar

23. Cheng S, Li Z, Gao R, Xing B, Gao Y, Yang Y, et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell. (2021) 184:792–809 e23. doi: 10.1016/j.cell.2021.01.010

PubMed Abstract | Crossref Full Text | Google Scholar

24. Krishna C, DiNatale RG, Kuo F, Srivastava RM, Vuong L, Chowell D, et al. Single-cell sequencing links multiregional immune landscapes and tissue-resident T cells in ccRCC to tumor topology and therapy efficacy. Cancer Cell. (2021) 39:662–77 e6. doi: 10.1016/j.ccell.2021.03.007

PubMed Abstract | Crossref Full Text | Google Scholar

25. Hong M, Tao S, Zhang L, Diao LT, Huang X, Huang S, et al. RNA sequencing: new technologies and applications in cancer research. J Hematol Oncol. (2020) 13:166. doi: 10.1186/s13045-020-01005-x

PubMed Abstract | Crossref Full Text | Google Scholar

26. Li X, Poire A, Jeong KJ, Zhang D, Chen G, Sun C, et al. Single-cell trajectory analysis reveals a CD9 positive state to contribute to exit from stem cell-like and embryonic diapause states and transit to drug-resistant states. Cell Death Discov. (2023) 9:285. doi: 10.1038/s41420-023-01586-9

PubMed Abstract | Crossref Full Text | Google Scholar

27. Smolander J, Junttila S, and Elo LL. Cell-connectivity-guided trajectory inference from single-cell data. Bioinformatics. (2023) 39:btad515. doi: 10.1093/bioinformatics/btad515

PubMed Abstract | Crossref Full Text | Google Scholar

28. Lu Y, Yang A, Quan C, Pan Y, Zhang H, Li Y, et al. A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma. Nat Commun. (2022) 13:4594. doi: 10.1038/s41467-022-32283-3

PubMed Abstract | Crossref Full Text | Google Scholar

29. Roessler S, Jia HL, Budhu A, Forgues M, Ye QH, Lee JS, et al. A unique metastasis gene signature enables prediction of tumor relapse in early-stage hepatocellular carcinoma patients. Cancer Res. (2010) 70:10202–12. doi: 10.1158/0008-5472.CAN-10-2607

PubMed Abstract | Crossref Full Text | Google Scholar

30. Grinchuk OV, Yenamandra SP, Iyer R, Singh M, Lee HK, Lim KH, et al. Tumor-adjacent tissue co-expression profile analysis reveals pro-oncogenic ribosomal gene signature for prognosis of resectable hepatocellular carcinoma. Mol Oncol. (2018) 12:89–113. doi: 10.1002/1878-0261.12153

PubMed Abstract | Crossref Full Text | Google Scholar

31. Jiang Z, Wu Y, Miao Y, Deng K, Yang F, Xu S, et al. HCCDB v2.0: decompose expression variations by single-cell RNA-seq and spatial transcriptomics in HCC. Genomics Proteomics Bioinf. (2024) 22:qzae011. doi: 10.1093/gpbjnl/qzae011

PubMed Abstract | Crossref Full Text | Google Scholar

32. Butler A, Hoffman P, Smibert P, Papalexi E, and Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. (2018) 36:411–20. doi: 10.1038/nbt.4096

PubMed Abstract | Crossref Full Text | Google Scholar

33. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. (2019) 16:1289–96. doi: 10.1038/s41592-019-0619-0

PubMed Abstract | Crossref Full Text | Google Scholar

34. Tran HTN, Ang KS, Chevrier M, Zhang X, Lee NYS, Goh M, et al. A benchmark of batch-effect correction methods for single-cell RNA sequencing data. Genome Biol. (2020) 21:12. doi: 10.1186/s13059-019-1850-9

PubMed Abstract | Crossref Full Text | Google Scholar

35. Meng Y, Ye F, Nie P, Zhao Q, An L, Wang W, et al. Immunosuppressive CD10(+)ALPL(+) neutrophils promote resistance to anti-PD-1 therapy in HCC by mediating irreversible exhaustion of T cells. J Hepatol. (2023) 79:1435–49. doi: 10.1016/j.jhep.2023.08.024

PubMed Abstract | Crossref Full Text | Google Scholar

36. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402

PubMed Abstract | Crossref Full Text | Google Scholar

37. Brunet JP, Tamayo P, Golub TR, and Mesirov JP. Metagenes and molecular pattern discovery using matrix factorization. Proc Natl Acad Sci U S A. (2004) 101:4164–9. doi: 10.1073/pnas.0308531101

PubMed Abstract | Crossref Full Text | Google Scholar

38. Hamamoto R, Takasawa K, Machino H, Kobayashi K, Takahashi S, Bolatkan A, et al. Application of non-negative matrix factorization in oncology: one approach for establishing precision medicine. Brief Bioinform. (2022) 23:bbac246. doi: 10.1093/bib/bbac246

PubMed Abstract | Crossref Full Text | Google Scholar

39. Lovmar L, Ahlford A, Jonsson M, and Syvanen AC. Silhouette scores for assessment of SNP genotype clusters. BMC Genomics. (2005) 6:35. doi: 10.1186/1471-2164-6-35

PubMed Abstract | Crossref Full Text | Google Scholar

40. Hanzelmann S, Castelo R, and Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | Crossref Full Text | Google Scholar

41. Zhang B and Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. (2005) 4:Article17. doi: 10.2202/1544-6115.1128

PubMed Abstract | Crossref Full Text | Google Scholar

42. Xu M, Zhou H, Hu P, Pan Y, Wang S, Liu L, et al. Identification and validation of immune and oxidative stress-related diagnostic markers for diabetic nephropathy by WGCNA and machine learning. Front Immunol. (2023) 14:1084531. doi: 10.3389/fimmu.2023.1084531

PubMed Abstract | Crossref Full Text | Google Scholar

43. Aran D, Hu Z, and Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | Crossref Full Text | Google Scholar

44. Mayakonda A, Lin DC, Assenov Y, Plass C, and Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | Crossref Full Text | Google Scholar

45. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | Crossref Full Text | Google Scholar

46. Maeser D, Gruener RF, and Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. (2021) 22:bbac246. doi: 10.1093/bib/bbab260

PubMed Abstract | Crossref Full Text | Google Scholar

47. Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X, et al. A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell. (2017) 171:1437–52 e17. doi: 10.1016/j.cell.2017.10.049

PubMed Abstract | Crossref Full Text | Google Scholar

48. Rebouissou S and Nault JC. Advances in molecular classification and precision oncology in hepatocellular carcinoma. J Hepatol. (2020) 72:215–29. doi: 10.1016/j.jhep.2019.08.017

PubMed Abstract | Crossref Full Text | Google Scholar

49. Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, et al. Integrated proteogenomic characterization of HBV-related hepatocellular carcinoma. Cell. (2019) 179:1240. doi: 10.1016/j.cell.2019.10.038

PubMed Abstract | Crossref Full Text | Google Scholar

50. Sharma A, Seow JJW, Dutertre CA, Pai R, Bleriot C, Mishra A, et al. Onco-fetal reprogramming of endothelial cells drives immunosuppressive macrophages in hepatocellular carcinoma. Cell. (2020) 183:377–94 e21. doi: 10.1016/j.cell.2020.08.040

PubMed Abstract | Crossref Full Text | Google Scholar

51. Jin H, Shi Y, Lv Y, Yuan S, Ramirez CFA, Lieftink C, et al. EGFR activation limits the response of liver cancer to lenvatinib. Nature. (2021) 595:730–4. doi: 10.1038/s41586-021-03741-7

PubMed Abstract | Crossref Full Text | Google Scholar

52. Rudalska R, Harbig J, Snaebjornsson MT, Klotz S, Zwirner S, Taranets L, et al. LXRalpha activation and Raf inhibition trigger lethal lipotoxicity in liver cancer. Nat Cancer. (2021) 2:201–17. doi: 10.1038/s43018-020-00168-3

PubMed Abstract | Crossref Full Text | Google Scholar

53. Shi Y, Zhang Q, Bi H, Lu M, Tan Y, Zou D, et al. Decoding the multicellular ecosystem of vena caval tumor thrombus in clear cell renal cell carcinoma by single-cell RNA sequencing. Genome Biol. (2022) 23:87. doi: 10.1186/s13059-022-02651-9

PubMed Abstract | Crossref Full Text | Google Scholar

54. Zhang K and Yuan E. Combined analysis of bulk and single-cell RNA sequencing reveals novel natural killer cell-related prognostic biomarkers for predicting immunotherapeutic response in hepatocellular carcinoma. Front Immunol. (2023) 14:1142126. doi: 10.3389/fimmu.2023.1142126

PubMed Abstract | Crossref Full Text | Google Scholar

55. Liu L, Liu Z, Gao J, Liu X, Weng S, Guo C, et al. CD8+ T cell trajectory subtypes decode tumor heterogeneity and provide treatment recommendations for hepatocellular carcinoma. Front Immunol. (2022) 13:964190. doi: 10.3389/fimmu.2022.964190

PubMed Abstract | Crossref Full Text | Google Scholar

56. Wang M, Han J, Xing H, Zhang H, Li Z, Liang L, et al. Dysregulated fatty acid metabolism in hepatocellular carcinoma. Hepat Oncol. (2016) 3:241–51. doi: 10.2217/hep-2016-0012

PubMed Abstract | Crossref Full Text | Google Scholar

57. Nomura M, Liu J, Rovira II, Gonzalez-Hurtado E, Lee J, Wolfgang MJ, et al. Fatty acid oxidation in macrophage polarization. Nat Immunol. (2016) 17:216–7. doi: 10.1038/ni.3366

PubMed Abstract | Crossref Full Text | Google Scholar

58. Su P, Wang Q, Bi E, Ma X, Liu L, Yang M, et al. Enhanced lipid accumulation and metabolism are required for the differentiation and activation of tumor-associated macrophages. Cancer Res. (2020) 80:1438–50. doi: 10.1158/0008-5472.CAN-19-2994

PubMed Abstract | Crossref Full Text | Google Scholar

59. Klebanoff CA, Gattinoni L, Torabi-Parizi P, Kerstann K, Cardones AR, Finkelstein SE, et al. Central memory self/tumor-reactive CD8+ T cells confer superior antitumor immunity compared with effector memory T cells. Proc Natl Acad Sci U S A. (2005) 102:9571–6. doi: 10.1073/pnas.0503726102

PubMed Abstract | Crossref Full Text | Google Scholar

60. Hunt EG, Hurst KE, Riesenberg BP, Kennedy AS, Gandy EJ, Andrews AM, et al. Acetyl-CoA carboxylase obstructs CD8(+) T cell lipid utilization in the tumor microenvironment. Cell Metab. (2024) 36:969–83 e10. doi: 10.1016/j.cmet.2024.02.009

PubMed Abstract | Crossref Full Text | Google Scholar

61. Scharping NE, Menk AV, Moreci RS, Whetstone RD, Dadey RE, Watkins SC, et al. The tumor microenvironment represses T cell mitochondrial biogenesis to drive intratumoral T cell metabolic insufficiency and dysfunction. Immunity. (2016) 45:374–88. doi: 10.1016/j.immuni.2016.07.009

PubMed Abstract | Crossref Full Text | Google Scholar

62. Chamoto K, Chowdhury PS, Kumar A, Sonomura K, Matsuda F, Fagarasan S, et al. Mitochondrial activation chemicals synergize with surface receptor PD-1 blockade for T cell-dependent antitumor activity. Proc Natl Acad Sci U S A. (2017) 114:E761–E70. doi: 10.1073/pnas.1620433114

PubMed Abstract | Crossref Full Text | Google Scholar

63. Liu X, Ren B, Ren J, Gu M, You L, and Zhao Y. The significant role of amino acid metabolic reprogramming in cancer. Cell Commun Signal. (2024) 22:380. doi: 10.1186/s12964-024-01760-1

PubMed Abstract | Crossref Full Text | Google Scholar

64. Shan X, Hu P, Ni L, Shen L, Zhang Y, Ji Z, et al. Serine metabolism orchestrates macrophage polarization by regulating the IGF1-p38 axis. Cell Mol Immunol. (2022) 19:1263–78. doi: 10.1038/s41423-022-00925-7

PubMed Abstract | Crossref Full Text | Google Scholar

65. Molotkov A, Deltour L, Foglio MH, Cuenca AE, and Duester G. Distinct retinoid metabolic functions for alcohol dehydrogenase genes Adh1 and Adh4 in protection against vitamin A toxicity or deficiency revealed in double null mutant mice. J Biol Chem. (2002) 277:13804–11. doi: 10.1074/jbc.M112039200

PubMed Abstract | Crossref Full Text | Google Scholar

66. Molotkov A, Fan X, and Duester G. Excessive vitamin A toxicity in mice genetically deficient in either alcohol dehydrogenase Adh1 or Adh3. Eur J Biochem. (2002) 269:2607–12. doi: 10.1046/j.1432-1033.2002.02935.x

PubMed Abstract | Crossref Full Text | Google Scholar

67. Long J, Wang A, Bai Y, Lin J, Yang X, Wang D, et al. Development and validation of a TP53-associated immune prognostic model for hepatocellular carcinoma. EBioMedicine. (2019) 42:363–74. doi: 10.1016/j.ebiom.2019.03.022

PubMed Abstract | Crossref Full Text | Google Scholar

68. Zhao H, Chen Q, Alam A, Cui J, Suen KC, Soo AP, et al. The role of osteopontin in the progression of solid organ tumour. Cell Death Dis. (2018) 9:356. doi: 10.1038/s41419-018-0391-6

PubMed Abstract | Crossref Full Text | Google Scholar

69. Tokuda K, Morine Y, Miyazaki K, Yamada S, Saito Y, Nishi M, et al. The interaction between cancer associated fibroblasts and tumor associated macrophages via the osteopontin pathway in the tumor microenvironment of hepatocellular carcinoma. Oncotarget. (2021) 12:333–43. doi: 10.18632/oncotarget.27881

PubMed Abstract | Crossref Full Text | Google Scholar

70. Nallasamy P, Nimmakayala RK, Karmakar S, Leon F, Seshachadryulu P, Lakshmanan I, et al. Pancreatic tumor microenvironment factor promotes cancer stemness via SPP1-CD44 axis. Gastroenterology. (2021) 161:1998–2013 e7. doi: 10.1053/j.gastro.2021.08.023

PubMed Abstract | Crossref Full Text | Google Scholar

71. Fako V, Martin SP, Pomyen Y, Budhu A, Chaisaingmongkol J, Franck S, et al. Gene signature predictive of hepatocellular carcinoma patient response to transarterial chemoembolization. Int J Biol Sci. (2019) 15:2654–63. doi: 10.7150/ijbs.39534

PubMed Abstract | Crossref Full Text | Google Scholar