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 |log2 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 (|log2FC| >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.
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 (IC50) 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).

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).

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:
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).

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.

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.

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).

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.

Supplementary
The article’s supplementary files as