Establishment of a prognostic model and immune profiling based on phagocytic regulatory genes in hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC) is the sixth most common cancer worldwide and accounts for the fourth highest number of cancer-related fatalities due to its significant mortality rate (1). Early-stage HCC can be effectively treated with surgical resection, liver transplantation, and local ablation (2,3). However, approximately 80% of HCC patients are diagnosed at an advanced stage, where the disease is either unresectable or has metastasized, leading to a poor prognosis (4). While targeted therapies combined with immunotherapy have showed some promise for advanced HCC patients, HCC remains one of the least responsive cancers to current immunotherapies (5). Therefore, further investigation into the HCC tumor microenvironment (TME) and the identification of effective biomarkers are urgently needed to enhance treatment efficacy for HCC patients.

HCC is an inflammation-associated tumor that is closely linked to immune tolerance and immune evasion within the TME (6). Tumor-associated macrophages (TAMs) represent the largest group of immune cells within the microenvironment of HCC, accounting for more than half of the immune cells present in the TME. They are crucial in fostering the immune suppression characteristic of HCC progression (7). TAMs are defined as macrophages that infiltrate tumor tissues or accumulate within the TME, and their polarization status is strongly associated with tumor initiation and progression (8,9). TAMs are classified into two distinct phenotypes: tumor-suppressive M1 macrophages and tumor-promoting M2 macrophages (8). Studies indicate that D-lactate has the potential to transform M2 TAMs into M1 TAMs through the PI3K/AKT signaling pathway, effectively improving survival in HCC mouse models (10-12). Consequently, TAMs have emerged as a promising therapeutic target in the treatment of HCC. It is worth mentioning that phagocytic regulatory factors (PRFs) play an important role in regulating the activity of tumor-induced immune cells such as macrophages and T cells (13,14). PRFs are expressed in nearly all cancer types and contribute to tumor development. One study demonstrated that the loss of IRF8 prevents the exhaustion of cancer-cell-reactive cytotoxic T lymphocytes and inhibits tumor growth (15). Overexpression of the surface receptor CD47 on TAMs stimulates macrophages in HCC to release immune evasion signals, thereby promoting tumor proliferation (16). These studies highlight the crucial role of TAM-associated PRFs in regulating tumor progression. However, there has been no systematic investigation on the prognostic value of phagocytic regulatory factor-related genes (PFRGs) in HCC. Therefore, a deeper understanding of the functional role of PFRGs in HCC would be highly beneficial in improving the efficacy of immunotherapy and the survival outcomes of HCC patients.

In this research, we divided The Cancer Genome Atlas (TCGA)-Liver Hepatocellular Carcinoma (LIHC) samples into two groups based on high and low PFRGs enrichment scores. Using this distinction, we then crafted a predictive model by pinpointing genes that vary in expression and are associated with HCC prognosis. We combined clinical data with these risk assessments to create a comprehensive nomogram. Moreover, we performed an in-depth analysis of the immune profile of HCC and pinpointed potential treatment options, providing valuable insights into the prognosis and potential therapeutic strategies for HCC patients. A complete list of abbreviations used throughout this manuscript is provided in Table S1. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-207/rc).

Methods

Data collection for analysis

The mRNA expression data of HCC-associated TCGA-LIHC samples (50 normal and 374 tumor samples) were obtained from TCGA database (https://portal.gdc.cancer.gov/). The latest data on copy number variations and relevant clinical information were also collected. Cases lacking survival data were excluded from this investigation. In addition, microarray data from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), specifically GSE76427, were downloaded to serve as a validation cohort. A total of 171 PFRGs were selected for further analysis based on insights gathered from published studies (17). This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Identification of differentially expressed genes related to PFRs

An analysis of differential gene expression between normal and tumor groups in the TCGA-LIHC dataset was conducted using the edgeR package, establishing criteria of |log₂ fold change (FC)| >0.585 and a false discovery rate (FDR) of less than 0.05. Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out on the differentially expressed genes (DEGs). The intersection of DEGs with the PFRGs was used to identify differentially expressed genes related to phagocytic regulatory factors (DEPFRGs).

Single-sample gene set enrichment analysis (ssGSEA)

To comprehensively investigate the potential biological mechanisms and biomarkers in TCGA-LIHC samples, we conducted ssGSEA on the DEPFRGs gene set using the GSVA package. The ssGSEA scores of TCGA-LIHC samples for the DEPFRGs gene set were calculated, and samples were classified into high and low score groups based on the median score. Kaplan-Meier (K-M) survival analysis was performed using the survival package. The estimate package was used to assess the Immune Score, Stromal Score, ESTIMATE Score, and Tumor Purity in the high and low score groups. Immune infiltration levels in both groups were quantified using the microenvironment populations (MCP)-counter and EPIC algorithms. In addition, we compiled and statistically analyzed the clinical characteristics of samples with complete clinical information from both groups. Differential analysis (|log₂FC| >1, FDR <0.05) and functional enrichment were performed between the high and low score groups.

The identification of prognostic biomarkers and the construction of a prognostic risk model for HCC

To further refine the prognostic model, we integrated clinical data with the expression levels of DEGs from high- and low-score groups, focusing on tumor samples from patients with survival time greater than 30 days. Univariate regression analysis of the genes was conducted using the survival package, and prognostic-related genes were identified based on a P value threshold of <0.01. To mitigate the risk of model overfitting, least absolute shrinkage and selection operator (LASSO) regression was applied to the candidate genes using the glmnet package, with cross-validation employed to select the optimal penalty parameter (lambda). This approach helped eliminate highly correlated genes, thus reducing model complexity. Subsequently, multivariate regression analysis was performed using the survival package.

Riskscore=∑βgene×Expgene[1]

Expgene indicates the expression levels of the selected genes and the βgene represents their corresponding risk coefficients.

The median risk score served as the threshold to divide the samples into high-risk and low-risk categories. Subsequently, survival curves were created based on these risk scores. To assess the predictive capability, we plotted receiver operating characteristic (ROC) curves utilizing the timeROC package, and computed the area under the curve (AUC) for survival at 1, 3, and 5 years. Additionally, we illustrated the distribution of risk scores and survival outcomes within the two groups. The model was validated using the GEO dataset. Finally, a K-M analysis was performed to produce survival curves for the identified feature genes.

Gene set enrichment analysis (GSEA)

Enrichment analysis was performed for the high- and low-risk groups using GSEA 4.3.2 software. GO and KEGG enrichment analyses were performed on the DEGs between the high- and low-risk groups for functional annotation and pathway analysis.

Develop and assess an integrated nomogram

Based on the risk group classification, the clinical characteristics of TCGA-LIHC patients were analyzed, and a baseline table was constructed using the flextable package. Univariate and multivariate regression analyses were conducted, incorporating both clinical data and prognostic model risk scores. A nomogram was developed using the rms package to predict patient survival rates at 1, 3, and 5 years. Calibration curves and decision curve analysis (DCA) were plotted to assess the predictive accuracy of the nomogram. A correlation analysis between clinical characteristics and prognostic risk scores was conducted for TCGA-LIHC patients. Patients were stratified into subgroups based on various clinical features, and K-M survival curves were subsequently plotted for each subgroup.

Immune landscape analysis and immune response prediction

We used the ssGSEA algorithm to gauge the infiltration levels of 29 distinct immune cell types. For each risk group, we calculated the scores with the ESTIMATE algorithm and then compared the results. The CIBERSORT algorithm was used to evaluate the immune infiltration status across various risk groups. We also analyzed the expression levels of immune checkpoints and generated boxplots to visualize the differences among risk groups. We obtained the immunophenoscore (IPS) scores for TCGA-LIHC from The Cancer Immunome Atlas (TCIA; https://tcia.at). Programmed cell death protein 1 (PD-1) and programmed cell death ligand 1 (PD-L1) checkpoint inhibitors have emerged as important targets in immunotherapy. To evaluate how well risk scores can predict responses to these treatments, we employed the IMvigor210CoreBiologies package to obtain transcriptomic and clinical data from patients in the IMvigor210 cohort treated with anti-PD-L1 therapy.

Tumor mutational burden (TMB) and antitumor drug analysis

Mutation data from TCGA-LIHC were used to analyze the mutations in the top 20 genes across risk groups. The results were visualized through a waterfall plot generated using the GenVisR package. To identify new therapeutic targets and enhance drug efficacy, we employed the CellMiner database (https://discover.nci.nih.gov/cellminer/) to assess antitumor drugs exhibiting a strong correlation with prognostic feature genes. To estimate the half-maximal inhibitory concentration (IC₅₀) of different medications across both high- and low-risk cohorts, we utilized the pRRophetic package. Notably, a lower IC50 value indicates greater drug efficacy against tumors.

Statistical analysis

All data analyses were performed using R software (version 4.4.1) and its associated packages. Statistical differences between the two groups were assessed by Wilcoxon test and t-test. K-M survival curves were analysed using the log-rank test to calculate P values between groups. The ROC curves were plotted using the timeROC package, and data visualisation relied heavily on the ggplot2 package. Correlations between variables were analysed using Spearman and Pearson methods. The median was used as the cut-off value for grouping. The statistical significance level was set at P<0.05.

Results

Identification and preliminary analysis of differentially expressed PFRGs

Initially, we performed differential expression analysis between the tumor and control groups in the TCGA-LIHC cohort, identifying 3,477 DEGs (Figure 1A). By intersecting the DEGs with the PFRGs gene set, we identified 19 DEPFRGs (Figure 1B). Notably, the majority of these DEPFRGs were found to be downregulated in HCC tumors (Figure 1C). To further explore the potential biological functions of these DEPFRGs, we conducted GO and KEGG enrichment analyses. The results indicated that these genes are primarily involved in biological processes such as transmembrane transport activity and extracellular matrix organization (Figure 1D). Furthermore, key pathways associated with these genes include cytokine receptor interactions and DNA replication (Figure 1E).

Figure 1 Identification and analysis of intersecting genes related to PFRGs. (A) A volcano plot highlights the upregulated and downregulated genes in the TCGA-LIHC cohort. (B) An upset plot demonstrates the overlap between DEGs and PFRGs, revealing 19 shared genes. The left bar represents the 3,458 DEGs remaining after the intersection analysis, the middle bar represents the 152 PFRGs remaining after the intersection analysis, and the right bar represents the 19 intersecting genes resulting from the intersection analysis performed on DEGs and PFRGs. (C) A boxplot illustrates the expression differences of the 19 DEPFRGs between normal and tumor samples. (D) GO enrichment analysis of DEPFRGs, and (E) KEGG enrichment analysis results. ***, P<0.001. DEG, differentially expressed gene; DEPFRGs, differentially expressed genes related to phagocytic regulatory factors; FC, fold change; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LIHC, Liver Hepatocellular Carcinoma; PFRGs, phagocytic regulatory factor-related genes; TCGA, The Cancer Genome Atlas.

Feature of PRFs enrichment in TCGA-LIHC

To gain a deeper understanding of the molecular characteristics and heterogeneity of TCGA-LIHC samples, we utilized the ssGSEA algorithm to assess the enrichment of the DEPFRGs gene set in these samples, generating enrichment scores associated with PRFs. Based on these scores, we stratified the samples into high- and low-scoring groups. The clinical baseline characteristics of HCC patients in these groups are presented in Table S2. Heatmaps and boxplots were created to visualize the differential expression of genes between the high- and low-scoring groups (Figure S1A and Figure 2A), with most genes showing higher expression levels in the high-scoring group. K-M survival analysis revealed significantly better survival in the high-scoring group compared to the low-scoring group (Figure 2B). Further examination of the various scores showed a notably greater tumor purity in the low-score group compared to the high-score group (Figure 2C). These findings suggest that PFRGs can effectively differentiate patients within the TCGA-LIHC cohort, highlighting their potential role in HCC progression. To further investigate immune cell infiltration, we employed the MCP-counter and EPIC algorithms, which revealed a higher level of T cell infiltration in the high-scoring group (Figure 2D,2E). Differential analysis between the two score groups identified 124 DEGs, with 106 upregulated and 18 downregulated genes (Figure S1B). Further investigation into the potential molecular mechanisms underlying the enrichment score groups revealed that these DEGs are critically involved in the regulation of T cell activation (Figure S1C). The primary pathway implicated is the cell adhesion molecule pathway (Figure S1D).

Figure 2 Feature of PRFs enrichment in TCGA-LIHC cohort. (A) The boxplot displays the expression differences of DEPFRGs between the two scoring groups. (B) K-M survival analysis comparing survival differences between the scoring groups. (C) Various scores across the scoring groups. (D) Assessment of immune cell infiltration levels in the scoring groups using the MCP-counter algorithm. (E) Evaluation of immune infiltration differences between the scoring groups using the EPIC algorithm. ns, P>0.5; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. CAF, consensus average function; DEPFRGs, differentially expressed genes related to phagocytic regulatory factors; EPIC, estimating the proportion of immune and cancer cells; K-M, Kaplan-Meier; LIHC, Liver Hepatocellular Carcinoma; MCP, microenvironment populations; NK, natural killer; PRFs, phagocytic regulatory factors; TCGA, The Cancer Genome Atlas.

Construction and validation of the prognostic risk model for HCC

A univariate regression analysis was conducted, leading to the identification of 26 genes potentially associated with HCC prognosis (Figure 3A, Table S3). To avoid model overfitting, LASSO analysis was employed, resulting in the selection of 9 feature genes (Figure S2A,S2B, Table S4). A multivariate regression analysis of these 9 genes was performed, and the final prognostic model was constructed based on them (Figure 3B, Table S5). The risk score was calculated as follows:

Riskscore=−0.0508∗CD5L−0.0178∗SLA2−0.0165∗FLT3−0.0749∗GPR18−0.0720∗BCL11B+0.0909∗CTSV−0.0113∗UBASH3A−0.0101∗XCR1−0.1351∗KLRK1[2]

Figure 3 Establishment and validation of a prognostic risk model for HCC. (A) Univariate regression analysis results used for identifying prognostic candidate genes. (B) Outcomes from multivariate regression analysis used for identifying prognostic feature genes. (C) Boxplot illustrating the expression differences of prognostic feature genes between the normal and tumor groups. (D) Boxplot showing the differential expression of prognostic feature genes between the risk groups. In the training set, the results from (E) ROC curves, (F) K-M curves, and (G) survival status distribution plots are presented. (H) In the validation set, the AUC values for 1, 3, and 5 years from the ROC analysis were all greater than 0.75. (I) The K-M curve indicates that the low-risk group demonstrates a better overall survival rate. (J) In GSE76427, risk score distribution and survival status distribution plots. ns, P>0.5; *, P<0.05; **, P<0.01; ***, P<0.001. AIC, akaike informat ion criterion; AUC, area under the curve; CI, confidence interval; K-M, Kaplan-Meier; ROC, receiver operating characteristic.

BCL11B, CD5L, FLT3, GPR18, KLRK1, SLA2, and XCR1 were significantly downregulated in the tumor group, whereas CTSV was found to be highly expressed there (Figure 3C). These genes were notably upregulated in the low-risk group (Figure 3D). To further assess the prognostic performance of the risk model, we calculated the AUC values for 1, 3, and 5 years using ROC curves, which were 0.701, 0.725, and 0.704, respectively (Figure 3E). Survival analysis indicated a better prognosis for the low-risk category (Figure 3F), with better survival outcomes compared to the high-risk group (Figure 3G). Additional validation with the GSE76427 gene set confirmed AUC values of 0.77, 0.787, and 0.886 for 1, 3, and 5 years, respectively (Figure 3H). K-M survival curves further supported these results, showing that the patients in the low-risk group had a higher overall survival rate and better survival outcomes (Figure 3I,3J). The findings of the K-M curve for the feature genes are shown in Figure S2C. These findings indicate that the risk model we constructed has robust prognostic performance.

Exploration of potential biological functions and active pathways in HCC

We conducted GSEA to further investigate the potential biological mechanisms in HCC. The high-risk group was predominantly enriched in biological processes related to DNA repair, RNA splicing, and cell differentiation (Figure 4A). In contrast, the low-risk group showed significant enrichment in processes associated with immune responses (Figure 4B). Furthermore, the DEGs between the risk groups were primarily associated with the regulation of immune cells (Figure 4C). The major pathways identified were linked to immune deficiencies and T cell activation (Figure 4D).

Figure 4 Investigation of the biological mechanisms and significant pathways in HCC. GSEA pathway enrichment analysis of the (A) high- and (B) low-risk groups. (C) GO and (D) KEGG enrichment analysis of DEGs between the high- and low-risk categories. DEGs, differentially expressed genes; FC, fold change; GSEA, gene set enrichment analysis; GO, Gene Ontology; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Development and evaluation of an integrated nomogram combining risk scores and clinical information for prognostic prediction in HCC patients

We conducted univariate regression analysis incorporating clinical factors and risk scores and identified that the prognostic factors, risk score (P<0.001), stage (P<0.001), and M stage (P=0.01), were statistically significant (Figure 5A). Figure 5B indicated that both risk score and overall stage serve as independent prognostic factors for HCC patients. Based on these findings, we included independent prognostic factors that were significant in both univariate and multivariate regression analyses: risk score and overall stage, and developed a comprehensive nomogram to enable researchers to efficiently evaluate the prognosis of HCC patients, thereby offering valuable support for clinical decision-making (Figure 5C). The accuracy of the nomogram was assessed through calibration curves, which demonstrated that the predicted overall survival rates at 1, 3, and 5 years closely aligned with the actual observed survival rates, underscoring the nomogram’s superior predictive performance (Figure 5D). Additionally, the DCA results highlighted that the nomogram provided enhanced clinical benefits across various threshold settings (Figure 5E). These findings suggested that our risk model could be well integrated with clinical features to help identify HCC patients in different risk tiers, improving the sensitivity and specificity of prognostic prediction to some extent. The model’s strong calibration and significant clinical decision-support capabilities enhanced its reliability in guiding treatment selection and prognostic management in real-world clinical settings.

Figure 5 Establishment and evaluation of a nomogram. (A) Univariate and (B) multivariate regression analysis based on clinical information and risk scores of HCC patients. (C) A nomogram with independent prognostic significance. (D) Calibration curve of the nomogram used to assess its accuracy. (E) The DCA curve showing the net benefit at various thresholds, providing valuable support for clinical decision-making. CI, confidence interval; DCA, decision curve analysis; HCC, hepatocellular carcinoma; M, metastasis; N, node.

Subgroup comparison analysis based on different clinical characteristics

Based on the risk assessment, HCC patients were categorized into high-risk and low-risk groups. The initial medical profiles for each group are depicted in Table 1. A comparison of risk scores between these subgroups revealed significant differences in risk scores across different T stages and pathological stages (Figure S3A). This suggested that the risk score we calculated effectively stratifies HCC patients based on their clinical features. Furthermore, we evaluated the overall survival of patients within these subgroups using K-M curves. Patients in the lower-risk group exhibited more favorable overall survival across all subgroups (Figure S3B).

Preliminary analysis of the immune landscape and immune therapy response

Immunotherapy has emerged as a promising approach in the treatment of HCC. In order to investigate the distribution of immune landscape across various risk levels, we utilized the ssGSEA algorithm to evaluate immune cell infiltration within two distinct risk groups. Our results revealed that most immune cells were markedly greater in the low-risk cohort (Figure S4). Additionally, the ESTIMATE score, immune score, and stromal score were substantially higher in the low-risk group, while tumor purity was notably lower in this cohort (Figure 6A). Moreover, the CIBERSORT algorithm showed significantly higher infiltration of M1 macrophages in the low-risk category (Figure 6B). M1 macrophages, recognized for their ability to suppress tumors, supported the hypothesis that individuals classified as low-risk are more likely to gain substantial advantages from immunotherapy. An examination of genes associated with immune checkpoints indicated considerably elevated expression levels within the low-risk cohort (Figure 6C). The IPS, a commonly used indicator for evaluating the efficacy of anti-PD-1 and CTLA-4 inhibitors, was notably elevated in the low-risk population, suggesting that these patients are more likely to exhibit a robust immune response to immunotherapy (Figure 6D). We evaluated the efficacy of PD-L1 blockade immunotherapy in the IMvigor210 cohort, which comprised 348 patients. Treatment responses included stable disease (SD), partial response (PR), complete response (CR), and disease progression (PD). Survival analysis revealed that low-risk patients had markedly improved overall survival compared to those in the high-risk group (Figure 6E). Individuals diagnosed with SD and PD were grouped as non-responders (NR), whereas those presenting with CR and PR were identified as responders (R). Interestingly, the response rate was markedly greater within the low-risk cohort (Figure 6F). Additionally, Figure 6G illustrated that responders demonstrated considerably lower risk scores compared to their non-responder counterparts. These findings suggested that patients with HCC who possess lower risk scores are more inclined to experience significant advantages from immunotherapy.

Figure 6 Assessment of the immune landscape and response to immunotherapy. (A) The ESTIMATE score, immune score, stromal score, and tumor purity metrics for various risk groups are presented. (B) The boxplot visualizes the disparities in the infiltration grades of several immune cells, applying the CIBERSORT approach. (C) The expression differences of immune checkpoint-related genes between risk groups. (D) The IPS in different risk stratified populations. (E) Overall survival analysis based on risk scores in the IMvigor210 cohort. (F) In the IMvigor210 cohort, the immune response status in high- and low-risk categories. (G) The cloud plot indicates the comparison of risk scores across different immune response groups. ns, P>0.05; *, P<0.05; ***, P<0.001. NR, non-responder group; IPS, immunophenoscore; R, responder group.

TMB and drug sensitivity profiling

To explore the genomic mutation landscape across risk groups, we analyzed mutation data from the TCGA-LIHC cohort and assessed the mutation frequency of the top 20 genes in both risk groups. These results showed that the risk groups had different mutational profiles (Figure S5). Additionally, the pRRophetic algorithm was employed to evaluate the IC50 values of four commonly used drugs in the treatment of HCC patients (Figure 7A,7B). Notably, the low-risk group exhibited greater sensitivity to Belinostat, Dabrafenib, and Sorafenib. Importantly, we identified several drugs whose sensitivity was significantly associated with the expression levels of specific genes. For instance, Ribavirin exhibited a significant positive correlation with BCL11B expression, while Afatinib showed a similar positive correlation with CTSV expression (Figure 7C). Conversely, some drugs, such as Sonidegib and Lenvatinib, displayed a negative correlation with the expression of GPR18 (Figure 7D).

Figure 7 Evaluation of TMB and drug sensitivity. (A) Violin plot illustrating three drugs with higher sensitivity in the low-risk category. (B) Violin plot demonstrating that Docetaxel exhibits greater sensitivity in the high-risk group. Drugs whose sensitivity is significantly correlated with the expression levels of specific feature genes, as predicted by the CellMiner database: (C) exhibiting a significant positive correlation and (D) exhibiting a significant negative correlation. IC50, half-maximal inhibitory concentration; TMB, tumor mutational burden.

Discussion

In our research, we divided the TCGA-LIHC samples into two categories, high and low, based on their PFRG enrichment scores. A differential gene analysis of these groups pinpointed nine pivotal genes with a significant link to the prognosis of HCC, including CD5L, SLA2, FLT3, GPR18, BCL11B, CTSV, UBASH3A, XCR1, and KLRK1. Using these genes, we crafted a robust risk assessment model. In order to explore the differences in prognosis between patients in different risk groups, we performed a comprehensive risk score assessment of patients with different clinical stages. By combining patients’ risk scores with their clinical characteristics, we also created a comprehensive nomogram that shows a robust net benefit at different thresholds, improving clinicians’ ability to make informed decisions. The model showed that patients in the low-risk group had significantly better prognostic outcomes. Further, a comprehensive prediction of immune cell infiltration, response to immunotherapy, and sensitivity to drugs was performed for different risk groups, and it was found that the low-risk group had higher levels of immune infiltration and better response to immunotherapy. In addition, drug sensitivity analyses provided different treatment strategies for patients in different risk groups.

Characteristic genes identified in this study are closely linked to tumor formation and progression. Among these, CD5L is a soluble glycoprotein that plays a role in modulating macrophage activity during the pathogenesis of various infections and inflammatory processes (18). One study has reported that CD5L promotes the polarization of M2 macrophages by regulating autophagy-mediated upregulation of ID3 (19). Another investigation suggests that in intrahepatic cholangiocarcinoma, CD5L-mediated crosstalk with CD8⁺ T cells leads to the upregulation of CTLA4 expression on CD8⁺ T cells, thereby impairing immune responses (20). These findings suggest that CD5L not only influences macrophage polarization but also undermines immune responses. Nevertheless, the potential molecular mechanisms underlying the role of CD5L in HCC remain to be fully elucidated. Additionally, XCR1, which is expressed in both murine and human conventional dendritic cells, has demonstrated a remarkable ability to inhibit viral replication (21). TIM-3 inhibits production of the chemokine CXCL9 by XCR1⁺ classical dendritic cells (cDC1), thereby limiting antitumor immunity in mammary carcinomas (22). These characteristic genes have exhibited notable regulatory functions in cancer, indicating their potential prognostic biomarkers for HCC.

Immunotherapy has significantly altered the landscape of cancer treatment, sparking a resurgence in tumor immunology research. Therapeutic strategies such as adoptive cell transfer (ACT) and immune checkpoint inhibitors (ICIs) have shown sustained clinical success, offering promising new avenues for patients (23). However, only a subset of cancer patients experiences significant benefits from these therapies (23). For instance, patients with HCC often exhibit limited responses to immunotherapy (24). Immune cell infiltration in the TME plays a pivotal role in tumor progression, thus influencing the clinical prognosis of cancer patients. In our study, we observed significantly higher expression of M1 macrophages in the low-risk group. The functional activity of M1 macrophages is predominantly reliant on aerobic glycolysis and plays a crucial role in enhancing tumor immunity (25). Hao et al. reported that in the context of HCC, the suppression of APOC1 expression facilitates the transition of M2 macrophages into M1 macrophages through a process known as ferroptosis. This shift contributes to the reconfiguration of the TME and enhances the effectiveness of anti-PD-1 immunotherapy (26). ICIs have now become the preferred first-line therapy for HCC (27), although numerous challenges remain. Our findings suggest that several immune checkpoint-related genes are significantly upregulated in the low-risk group, indicating that these patients may derive greater benefits from immunotherapy. Moreover, one study has shown that in HCC, CCL19 recruits dendritic cells (DCs) to interact with CD4⁺ T cells, thereby promoting the differentiation of CD4⁺ T cell subtypes and enhancing the anti-tumor response (28). HCC patients with reduced expression of TNFSF14 typically exhibit more aggressive tumors and poorer prognosis. The loss of TNFSF14 impairs the apoptosis of HepG2 cells induced by LX-2 and NaHS, a process mediated by the JNK/JunB signaling pathway (29). Functional enrichment analysis conducted in this study revealed that the DEGs between high- and low-risk groups are primarily involved in immune regulation pathways. These findings underscore the significant role of immune modulation in the treatment of HCC. Moreover, the extent of immune cell infiltration and the expression levels of immune checkpoint-related genes offer novel perspectives for reshaping the HCC immune microenvironment and refining immunotherapeutic strategies for HCC.

In this study, we identified three anti-tumor drugs, belinostat, dabrafenib, and sorafenib, that exhibit greater sensitivity in the low-risk cohort. Belinostat is an effective histone deacetylase inhibitor (HDAC) that induces cytotoxicity in Kirsten rat sarcoma virus-induced lung cancer cells by modulating metabolic reprogramming (30). Dabrafenib, an ATP-competitive Raf inhibitor, and sorafenib, an orally active Raf inhibitor, both function as multi-kinase inhibitors, inducing autophagy and apoptosis in tumor cells, thus demonstrating anti-tumor effects (31,32). Sorafenib has been investigated in the context of HCC treatment, with research by Lu et al. showing that Sorafenib resistance in HCC is mediated by the ETS1/miR-23a-3p/ACSL4 axis, which regulates ferroptosis (33). Additionally, we found that Docetaxel exhibited greater sensitivity in the high-risk group. Docetaxel is a microtubule depolymerization agent that induces apoptosis and can remodel the immune microenvironment in prostate cancer, enhancing ICI-based therapies (34). The broad clinical application of these drugs in cancer treatment suggests their potential as therapeutic agents for HCC, with significant clinical value. Notably, we identified four drugs, ribavirin, afatinib, sonidegib, and lenvatinib, that are significantly associated with the expression of feature genes. These drugs exhibit potential therapeutic value in the treatment of HCC and merit further investigation.

In summary, our study investigates the PRFs associated with TAMs in HCC, offering new insights into their crucial roles in HCC prognosis and the immune microenvironment. By constructing a prognostic risk assessment model based on nine feature genes, we provide valuable references and recommendations for risk stratification and clinical management of HCC patients.

Conclusions

In conclusion, this study identified nine prognostic biomarkers associated with HCC by analyzing PFRGs, and constructed a risk assessment model. The findings suggest that HCC patients with lower risk scores have a better prognosis. The nomogram, developed by integrating clinical information and risk scores, offers important insights for clinical decision-making. Additionally, we systematically assessed the immune landscape of HCC, offering new perspectives for immunotherapy in HCC.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-207/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-207/prf

Funding: The study was approved by 2022 Tangshan Municipal Science and Technology Plan Project (22130222G).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-207/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Brown ZJ, Tsilimigras DI, Ruff SM, et al. Management of Hepatocellular Carcinoma: A Review. JAMA Surg 2023;158:410-20. [Crossref] [PubMed]
2. Xiao Y, Chen J, Wang J, et al. Acute Myeloid Leukemia Epigenetic Immune Escape From Nature Killer Cells by ICAM-1. Front Oncol 2021;11:751834. [Crossref] [PubMed]
3. Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA Cancer J Clin 2024;74:12-49. [Crossref] [PubMed]
4. Vogel A, Meyer T, Sapisochin G, et al. Hepatocellular carcinoma. Lancet 2022;400:1345-62. [Crossref] [PubMed]
5. Zeng L, Su J, Qiu W, et al. Survival Outcomes and Safety of Programmed Cell Death/Programmed Cell Death Ligand 1 Inhibitors for Unresectable Hepatocellular Carcinoma: Result From Phase III Trials. Cancer Control 2022;29:10732748221092924. [Crossref] [PubMed]
6. Zhang Y, Han G, Gu J, et al. Role of tumor-associated macrophages in hepatocellular carcinoma: impact, mechanism, and therapy. Front Immunol 2024;15:1429812. [Crossref] [PubMed]
7. Huang Y, Ge W, Zhou J, et al. The Role of Tumor Associated Macrophages in Hepatocellular Carcinoma. J Cancer 2021;12:1284-94. [Crossref] [PubMed]
8. Cheng K, Cai N, Zhu J, et al. Tumor-associated macrophages in liver cancer: From mechanisms to therapy. Cancer Commun (Lond) 2022;42:1112-40. [Crossref] [PubMed]
9. Boutilier AJ, Elsawa SF. Macrophage Polarization States in the Tumor Microenvironment. Int J Mol Sci 2021;22:6995. [Crossref] [PubMed]
10. Han S, Bao X, Zou Y, et al. d-lactate modulates M2 tumor-associated macrophages and remodels immunosuppressive tumor microenvironment for hepatocellular carcinoma. Sci Adv 2023;9:eadg2697. [Crossref] [PubMed]
11. McDonald B, Zucoloto AZ, Yu IL, et al. Programing of an Intravascular Immune Firewall by the Gut Microbiota Protects against Pathogen Dissemination during Infection. Cell Host Microbe 2020;28:660-668.e4. [Crossref] [PubMed]
12. Fruman DA, Chiu H, Hopkins BD, et al. The PI3K Pathway in Human Disease. Cell 2017;170:605-35. [Crossref] [PubMed]
13. Zhang W, Zhang Q, Yang N, et al. Crosstalk between IL-15Rα(+) tumor-associated macrophages and breast cancer cells reduces CD8(+) T cell recruitment. Cancer Commun (Lond) 2022;42:536-57. [Crossref] [PubMed]
14. Chen S, Wang M, Lu T, et al. JMJD6 in tumor-associated macrophage regulates macrophage polarization and cancer progression via STAT3/IL-10 axis. Oncogene 2023;42:2737-50. [Crossref] [PubMed]
15. Nixon BG, Kuo F, Ji L, et al. Tumor-associated macrophages expressing the transcription factor IRF8 promote T cell exhaustion in cancer. Immunity 2022;55:2044-2058.e5. [Crossref] [PubMed]
16. Pan Y, Yin Q, Wang Z, et al. AFP shields hepatocellular carcinoma from macrophage phagocytosis by regulating HuR-mediated CD47 translocation in cellular membrane. Transl Oncol 2025;52:102240. [Crossref] [PubMed]
17. Xin S, Sun X, Jin L, et al. The Prognostic Signature and Therapeutic Value of Phagocytic Regulatory Factors in Prostate Adenocarcinoma (PRAD). Front Genet 2022;13:877278. [Crossref] [PubMed]
18. Sanchez-Moral L, Paul T, Martori C, et al. Macrophage CD5L is a target for cancer immunotherapy. EBioMedicine 2023;91:104555. [Crossref] [PubMed]
19. Watanabe Y, Fukuda T, Hayashi C, et al. Extracellular vesicles derived from GMSCs stimulated with TNF-α and IFN-α promote M2 macrophage polarization via enhanced CD73 and CD5L expression. Sci Rep 2022;12:13344. [Crossref] [PubMed]
20. Lu JC, Wu LL, Sun YN, et al. Macro CD5L(+) deteriorates CD8(+)T cells exhaustion and impairs combination of Gemcitabine-Oxaliplatin-Lenvatinib-anti-PD1 therapy in intrahepatic cholangiocarcinoma. Nat Commun 2024;15:621. [Crossref] [PubMed]
21. Heger L, Hatscher L, Liang C, et al. XCR1 expression distinguishes human conventional dendritic cell type 1 with full effector functions from their immediate precursors. Proc Natl Acad Sci U S A 2023;120:e2300343120. [Crossref] [PubMed]
22. de Mingo Pulido Á, Hänggi K, Celias DP, et al. The inhibitory receptor TIM-3 limits activation of the cGAS-STING pathway in intra-tumoral dendritic cells by suppressing extracellular DNA uptake. Immunity 2021;54:1154-1167.e7. [Crossref] [PubMed]
23. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol 2020;17:807-21. [Crossref] [PubMed]
24. Sangro B, Sarobe P, Hervás-Stubbs S, et al. Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2021;18:525-43. [Crossref] [PubMed]
25. Li M, Yang Y, Xiong L, et al. Metabolism, metabolites, and macrophages in cancer. J Hematol Oncol 2023;16:80. [Crossref] [PubMed]
26. Hao X, Zheng Z, Liu H, et al. Inhibition of APOC1 promotes the transformation of M2 into M1 macrophages via the ferroptosis pathway and enhances anti-PD1 immunotherapy in hepatocellular carcinoma based on single-cell RNA sequencing. Redox Biol 2022;56:102463. [Crossref] [PubMed]
27. Pinter M, Scheiner B, Pinato DJ. Immune checkpoint inhibitors in hepatocellular carcinoma: emerging challenges in clinical practice. Lancet Gastroenterol Hepatol 2023;8:760-70. [Crossref] [PubMed]
28. Lu LL, Xiao SX, Lin ZY, et al. GPC3-IL7-CCL19-CAR-T primes immune microenvironment reconstitution for hepatocellular carcinoma therapy. Cell Biol Toxicol 2023;39:3101-19. [Crossref] [PubMed]
29. Cheng J, Li Y, Wang X, et al. Response stratification in the first-line combined immunotherapy of hepatocellular carcinoma at genomic, transcriptional and immune repertoire levels. J Hepatocell Carcinoma 2021;8:1281-1295. [Crossref] [PubMed]
30. Peter RM, Sarwar MS, Mostafa SZ, et al. Histone deacetylase inhibitor belinostat regulates metabolic reprogramming in killing KRAS-mutant human lung cancer cells. Mol Carcinog 2023;62:1136-46. [Crossref] [PubMed]
31. Takano K, Munehira Y, Hatanaka M, et al. Discovery of a Novel ATP-Competitive MEK Inhibitor DS03090629 that Overcomes Resistance Conferred by BRAF Overexpression in BRAF-Mutated Melanoma. Mol Cancer Ther 2023;22:317-32. [Crossref] [PubMed]
32. El-Ashmawy NE, Khedr EG, El-Bahrawy HA, et al. Sorafenib effect on liver neoplastic changes in rats: more than a kinase inhibitor. Clin Exp Med 2017;17:185-91. [Crossref] [PubMed]
33. Lu Y, Chan YT, Tan HY, et al. Epigenetic regulation of ferroptosis via ETS1/miR-23a-3p/ACSL4 axis mediates sorafenib resistance in human hepatocellular carcinoma. J Exp Clin Cancer Res 2022;41:3. [Crossref] [PubMed]
34. Ma Z, Zhang W, Dong B, et al. Docetaxel remodels prostate cancer immune microenvironment and enhances checkpoint inhibitor-based immunotherapy. Theranostics 2022;12:4965-79. [Crossref] [PubMed]