Comprehensive bioinformatic analysis reveals sorafenib response-related prognostic signature in hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC), accounting for approximately 90% of liver cancers, is one of the most deadly malignant tumors with high morbidity and mortality. HCC is reported to be the second leading cause of cancer-related death worldwide, with a 5-year survival rate of about 10% (1). The treatments for HCC include radiofrequency ablation, resection, or liver transplantation, but these strategies are only effective when diagnosed early (2). Nevertheless, most patients are relatively asymptomatic at an early stage, and 80% patients are diagnosed at an advanced stage when treatment options are very limited. HCC is known to be resistant to chemotherapy (3), thus, it is urgent to search for alternative treatment strategies.

The molecular targeted drugs have been designed to antagonize “oncogenic operable genes” that are specifically altered in tumors (4). Although many molecular-targeted agents have been tested in HCC, only lenvatinib and sorafenib have been approved as first-line agents for advanced HCC (5). However, the lack of second-line treatment options after lenvatinib means that sorafenib is the recommended gold standard for the treatment of advanced HCC (6,7). Sorafenib is a kind of multi-kinase inhibitor that can inhibit the proliferation of tumor cells by inhibiting the activity of serine/threonine kinases and receptor tyrosine kinases on tumor cells and tumor blood vessels (8). Sorafenib can prolong survival in HCC patients, but its efficacy is poor due to the progression of resistant cells (9). Unfortunately, the mechanisms underlying individual differences in response to sorafenib have not been fully elucidated (10). Currently, some studies have been applied to identify markers predicting sorafenib resistance in HCC patients. For instance, FGF3/4 (11) and VEGFA (12) have been reported to be more frequently amplified in sorafenib responders. Yuan et al. (2) identified four sorafenib resistance-related genes (IGF2R, C2orf27A, PON1, and CFB) in HCC. However, due to the heterogeneity of HCC, these findings are far from sufficient to explain the molecular mechanisms of sorafenib resistance. Development of appropriate biomarkers for patient stratification of sorafenib response may help physicians guide the selection of tailored therapies.

In recent years, the role of tumor microenvironment (TME) in human cancers has attracted more and more attention (13,14). However, to date, there has been a lack of discussion regarding the cellular characteristics of immune infiltration and its relationship with sorafenib response in HCC patients treated with sorafenib. As a result, we intended to identify gene signatures associated with response to sorafenib in HCC and constructed a prognostic risk model based on these gene signatures. Furthermore, the correlation of risk grouping with immune cells was explored. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-1020/rc).

Methods

Data sources

The gene expression level data of HCC from The Cancer Genome Atlas (TCGA) were downloaded from Xena database (https://xena.ucsc.edu/), and the detection platform was Illumina HiSeq 2000 RNA Sequencing. A total of 423 samples were included, including 367 tumor samples with clinical prognosis information and 50 normal control samples. This dataset was used as the training dataset.

The following two datasets were also downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) dataset (https://www.ncbi.nlm.nih.gov/geo/):

GSE109211 (15): a total of 140 HCC samples were included, and 67 of them receiving sorafenib treatment (46 non-responders and 21 responders) were included in the analysis. The detection platform was GPL13938 Illumina HumanHT-12 WG-DASL V4.0 expression beadchip. This dataset was used as the data for the screening analysis of genes related to sorafenib treatment response.

GSE76427 (16): 167 HCC-related genes were included, of which 115 HCC tumor samples were selected for inclusion. The assay platform was GPL10558 Illumina HumanHT-12 V4.0 expression beadchip. This dataset was used as a validation dataset for the prognostic model. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

HCC sensitivity to sorafenib evaluation and its difference and prognostic correlation analysis

Based on the genome-wide expression level data of all samples included in TCGA analysis, combined with the drug information in the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) (17), R4.3.1 pRRophetic (18) was used to estimate sorafenib half maximal inhibitory concentration (IC₅₀) in each sample. Then, Kruskal-Wallis was used to analyze the difference of IC₅₀ levels between HCC and normal control samples. In addition, combined with the clinical prognostic information of HCC tumor samples, samples were divided into high and low IC₅₀ groups according to the median value of IC₅₀. The Kaplan-Meier curve method in R4.3.1 survival package version 2.41-1 (19) was used to evaluate the correlation between high/low IC₅₀ grouping and clinical survival prognosis.

Screening of modules associated with HCC sensitivity to sorafenib by Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA was conducted on the basis of the genome-wide expression level data of all samples in TCGA using WGCNA package 1.61 (20) in R4.3.1. The module partition threshold was: cutHeight =0.995. After the modules were obtained, the correlation between the IC₅₀ of HCC samples and each module was calculated. The modules with a significant positive correlation with the IC₅₀ level (the correlation coefficient higher than 0.3) were retained and the genes contained in the modules were as the genes related to HCC sensitivity to sorafenib.

Screening of genes associated with sorafenib response

TCGA and GSE109211 datasets were used in this step. In the TCGA, the differentially expressed genes (DEGs) in HCC vs. normal were selected under the thresholds of false discovery rate (FDR) <0.05 and |log₂fold change (FC)| >1 using Limma (21). In GSE109211, the DEGs between responder and non-responder were selected under the same thresholds. The overlapping genes were retained as the significant DEGs associated with sorafenib response in HCC. These genes were then compared with the genes associated with HCC sensitivity to sorafenib screened by WGCNA. The overlapping genes of the two parts were integrated as the genes associated with sorafenib response.

Construction of prognostic model based on sorafenib response-related gene

Based on the DEGs associated with sorafenib response screened above, combined with the clinical survival and prognosis information of tumor samples in TCGA dataset, the genes significantly associated with prognosis were screened by univariate Cox regression analysis in survival package 2.41-1 in R4.3.1. A P value less than 0.05 was selected as the threshold for screening significant associations.

The selected prognosis-related genes were performed least absolute shrinkage and selection operator (LASSO) regression analysis using lars package 1.2 (22). The optimal combination of genes associated with sorafenib response was screened.

Then, according to the LASSO prognostic coefficient of each element in the optimal related gene combination and the expression level of the target gene in the dataset, the following risk score (RS) model was constructed:

Riskscore=esum(eachgene'sexpressionlevels∗correspondingcoefficient)esum(eachgene'smeanexpressionlevels∗correspondingcoefficient)[1]

Evaluation of RS model

In TCGA training set and GSE76427 validation dataset, the RS values were calculated respectively, and then the samples were divided into high (RS score higher than or equal to the median RS value) and low (RS score lower than the median RS value) sample groups according to the median RS value as the boundary. The survival 2.41-1 was used to evaluate the correlation between the high/low grouping and the actual survival prognosis.

Construction of a nomogram model

In the TCGA training set tumor samples, independent prognostic clinical factors were screened through univariate and multivariate Cox regression analysis in survival 2.41-1 by combining with the clinical information of the samples. A log-rank P value of less than 0.05 was selected as the threshold.

To further investigate the prognostic independence between independent prognostic clinical factors and RS factors, the nomogram prediction model for 1-, 3- and 5-year survival rates was constructed using the rms package version 5.1-2 (23,24) based on the independent prognostic factors, combined with the risk information identified by the prognostic prediction model. After that, the C-index coefficient of the nomogram model was calculated by R4.3.1 survcomp 1.34.0 (25). In addition, decision curve analyses were performed separately for each independent prognostic factor using rmda 1.6.

Immune distribution assessment and comparison

CIBERSORT (26) was used to evaluate the proportion of immune cell types in TCGA tumor samples, and then the distribution of different immune cells in different risk groups was compared. R4.3.1 estimate package (27) was used to calculate the ESTIMATE score, immune score, stromal score and tumor purity of TCGA tumor samples, and the distribution difference of each score in different risk groups was also compared. Finally, the correlation between the expression levels of the model genes and the immune cells/evaluation scores with significant differences was calculated.

Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway related to prognostic risk grouping

Based on all genomic expression levels detected in TCGA samples, Gene Set Enrichment Analysis (GSEA) (28) was used to screen KEGG signaling pathways significantly related to risk grouping. In the GSEA analysis results, there were three key statistical values: enrichment score (ES); normalized enrichment score (NES); nominal P value. Combined with the above three key statistic values, we finally selected an FDR value less than 0.05 as the threshold for screening significantly enriched related KEGG pathways.

The analysis flow is shown in Figure S1.

Statistical analysis

All analysis was performed by packages in R4.3.1. The thresholds of FDR <0.05 or P value <0.05 was considered the level of statistical significance.

Results

HCC sensitivity to sorafenib evaluation and its difference and prognostic correlation analysis

The results showed that HCC tumor samples had a lower IC₅₀ value, indicating that HCC tumor samples had a higher sensitivity to sorafenib than normal control samples (Figure 1A). According to the median value of sorafenib IC₅₀, the HCC samples were divided into high and low sorafenib IC₅₀ groups. Kaplan-Meier (KM) curve showed that HCC samples with low IC₅₀ (high IC₅₀ sensitivity) had better survival prognosis (Figure 1B).

Figure 1 Sorafenib IC₅₀ evaluation and its difference and prognostic correlation analysis. (A) Distribution of IC₅₀ levels in HCC and control samples. (B) KM curves of survival prognosis grouped by different IC₅₀ levels. ***, P<0.001. CI, confidence interval; CTRL, control; HCC, hepatocellular carcinoma; IC₅₀, half maximal inhibitory concentration; KM, Kaplan-Meier.

Screening of modules associated with HCC sensitivity to sorafenib

The power value (power =20) when the square value of the correlation coefficient reached 0.9 for the first time was selected (Figure 2A). The average node connectivity degree of the co-expression network constructed was 1, which was completely consistent with the properties of the small-world network. Then, the dissimilarity coefficient between nodes was calculated, and the system clustering tree was obtained. The pruning height was set as cutHeight =0.995, and 7 modules were obtained (Figure 2B). Finally, the correlation between the sorafenib value of the sample and each module was calculated (Figure 2C). The yellow, blue, green, and brown modules were significantly associated with the HCC sensitivity to sorafenib of the samples, which contained 798 genes.

Figure 2 Screening of modules associated with HCC sensitivity to sorafenib. (A) The parameters of the network construction. Left: power selection diagram of adjacency matrix weight parameter. The horizontal axis represents the weight parameter power, and the vertical axis represents the squared correlation coefficient between log(k) and log(p(k)) in the corresponding network. The red line indicates the standard line where the squared value of the correlation coefficient reaches 0.9. Right: schematic diagram of the average connectivity degree of genes under different power parameters. The red line represents the value of the average connectivity degree of network nodes under the value of the power parameter of the adjacency matrix weight in the left figure. (B) Module partition tree diagram. Each color represents a different module. (C) Correlation heatmap between each module and HCC sensitivity to sorafenib. HCC, hepatocellular carcinoma; IC₅₀, half maximal inhibitory concentration.

Screening of genes associated with sorafenib response

In TCGA dataset, 1,297 DEGs were selected in HCC vs. normal groups; in GSE109211, 2,315 DEGs were identified between responder and non-responder groups. A total of 399 overlapping significantly DEGs were obtained (Figure 3A).

Figure 3 Screening of genes associated with sorafenib response. (A) Venn diagram comparing significantly differentially expressed genes in the two datasets. (B,C) Common differentially expressed genes associated with (B) GO BP and (C) KEGG pathways. (D) Venn diagram of significantly differentially expressed genes overlapped between the two data sets compared with genes in the WGCNA module. TCGA, The Cancer Genome Atlas; BP, biological process; DEGs, differentially expressed genes; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; WGCNA, Weighted Gene Co-Expression Network Analysis.

After that, functional analysis of the overlapping genes was carried out for the enrichment of Gene Ontology (GO) biological processes and KEGG pathways. A total of 16 significantly related biological processes (negative regulation of cell proliferation, positive regulation of apoptotic process, response to xenobiotic stimulus, etc.) and 14 KEGG pathways (metabolic pathway, complement and coagulation cascades, alcoholic liver disease, drug metabolism-cytochrome P450, etc.) were screened (Figure 3B,3C).

Finally, 399 overlapping DEGs were compared with 798 module genes screened by WGCNA. As shown in Figure 3D, 56 overlapping genes were screened, which were considered as sorafenib response-related genes.

Construction of prognostic model based on sorafenib response-related gene

On the basis of 56 sorafenib response-related genes, and the sample clinical prognostic information, 14 survival prognosis-related genes were obtained, as presented in forest plot (Figure 4A). After LASSO regression analysis, 5 optimal genes were obtained: deoxyribonuclease 1 like 3 (DNASE1L3), long-chain acyl-CoA synthetase 6 (ACSL6), aggrecan (ACAN), brain-specific kinases 1 (BRSK1) and CD68 (Figure 4B). The five genes were grouped according to their median expression levels, and the correlation between different expression levels of each gene and survival prognosis was analyzed by KM curve. The results are shown in Figure 4C.

Figure 4 Construction of prognostic model based on sorafenib response-related gene. (A) Univariate cox regression analysis forest plots. (B) LASSO parametric map. (C) KM curves for five sorafenib-associated DEGs, with blue and red indicating low and high expression groups, respectively. CI, confidence interval; DEGs, differentially expressed genes; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; KM, Kaplan-Meier.

Evaluation of RS model

In TCGA and GSE76427 datasets, the RS values were respectively calculated, and then the samples in two datasets were divided into high- and low-risk groups. Then, the KM curve was used to evaluate the association between the high- and low-risk grouping and the actual prognostic information of the disease. The KM curve, the RS distribution and survival time for individual, and the receiver operating characteristic (ROC) curve results based on the five optimal genes of the training set are shown in Figure 5A-5C, and those of the validation set GSE76427 are shown in Figure 5D-5F. The area under curve (AUC) values for 1-, 3- and 5-year were 0.828, 0.775 and 0.759 respectively in TCGA dataset. In GSE76427 dataset, the AUC values were also larger than 0.7. The results showed that in TCGA training set and GSE76427 dataset, there was a significant correlation between the different risk groups of samples based on RS model prediction and the actual prognosis.

Figure 5 Evaluation of RS model. (A,D) KM curve plots associated with the five signature gene prognostic model in TCGA training set and validation set. (B,E) RS distribution, survival time status in TCGA training set and validation set. (C,F) ROC curves based on prognostic model. AUC, area under curve; CI, confidence interval; HR, hazard ratio; KM, Kaplan-Meier; ROC, receiver operating characteristic; RS, risk score; TCGA, The Cancer Genome Atlas.

Screening of clinical independent prognostic factors

In the clinical information of tumor samples in TCGA training set, univariate Cox regression analysis was used to screen the factors significantly related to survival prognosis, and age, pathologic stage and RS model status were obtained. The forest display diagram is shown in Figure 6A. Then multivariate cox regression analysis was used to screen the clinical prognostic factors related to independent prognosis, and 2 independent prognostic factors were screened: pathologic stage and RS model status. The KM curves associated with pathologic stage are shown in Figure 6B.

Figure 6 Screening of clinical independent prognostic factors and construction of nomogram model. (A) Forest plot presentation of univariate cox regression analysis. (B) KM curves associated with pathologic stage. (C) Nomogram of independent prognostic factors for survival prediction model. (D) Consistency line of 1-, 3-, and 5-year survival prediction and actual survival rate (horizontal axis represents the predicted survival rate; vertical axis represents the actual survival rate). (E) Decision line analysis curve of the sample. CI, confidence interval; HR, hazard ratio; KM, Kaplan-Meier; M, metastasis; N, node; OS, overall survival; RS, risk score; T, tumor.

Construction of a nomogram model

To further analyze the correlation between pathologic stage/RS model status factors and survival prognosis, a nomogram survival model was constructed (Figure 6C). The nomogram integrated each clinical indicator through the “Total points” axis in the first row to predict the survival of the sample. The 1-, 3- and 5-year survival rates predicted by the nomogram survival model were consistent with the actual 1-, 3- and 5-year survival rates (Figure 6D). In addition, we also performed a decision line analysis of the individual factors, as shown in Figure 6E.

Correlation analysis of risk grouping and immunity

Using CIBERSORT and Kruskal-Wallis test, nine immune cell types were identified to be significantly differentially distributed between risk groups, including B cell naïve, T cell CD4⁺ memory resting, T cell follicular helper, T cell regulatory (Tregs), NK cell resting, monocyte, macrophage M0, mast cell activated, neutrophil (Figure 7A). Additionally, the differences of the estimate scores among different risk groups were significant (Figure 7B). The correlations between the five model genes and nine immune cells with different distribution/each score were calculated (Figure 7C). Furthermore, HCC sensitivity to sorafenib was negatively correlated with RS value, and the low-risk group had a higher HCC sensitivity to sorafenib (Figure 7D,7E).

Figure 7 Correlation analysis of risk grouping and immunity. (A) Distribution of various immune cell types in samples from different risk groups. (B) Distribution of each score in samples with different risk groups. ***, P<0.001. (C) Correlation scores of five genes with immune cells. (D) Distribution of sorafenib IC₅₀ levels in different risk groups is shown. (E) Scatter plot of the association between sorafenib IC₅₀ levels and RS risk values. IC₅₀, half maximal inhibitory concentration; NK, natural killer; RS, risk score.

KEGG pathway analysis

Based on the genome-wide expression levels in TCGA samples, GSEA was used to screen KEGG signaling pathways significantly related to subtype grouping. A total of 30 significant pathways were obtained, such as antigen processing and antigen presentation, ascorbate and aldarate metabolism, beta alanine metabolism, and NOD-like receptor signaling pathway (Table S1).

Discussion

In early-stage HCC, local ablation and surgical resection are the mainstay of treatment, while for patients with advanced unresectable HCC, sorafenib is the first-line treatment (29). However, sorafenib is effective in very few patients and is associated with serious adverse reactions and drug resistance. Therefore, there is an urgent need for biomarkers to determine which patients are more likely to benefit from sorafenib. In this study, 399 significantly DEGs associated with sorafenib response were obtained, which were involved in metabolic and complement pathways. After LASSO analysis, five signature genes (DNASE1L3, ACSL6, ACAN, BRSK1 and CD68) were identified, and a RS prognostic risk prediction model was constructed. The ROC curves suggested that the model presented high predictive precision. Additionally, a nomogram was constructed based on pathologic stage and RS model status, which also had good prediction performance in survival prognosis. Moreover, we found significant correlation between risk groups and immunity. The patients in the high risk group may be benefit from sorafenib treatment.

Sorafenib, a multi-kinase inhibitor, has been the most widely used agent for the treatment of HCC patients at an advanced stage since 2007 (30). Our study found that HCC patients with lower sorafenib IC₅₀ value (higher sensitivity) showed better survival prognosis compared with normal controls, which directly displayed the positive correlation between responsiveness and clinical prognosis. Further analysis showed that patients in the high-risk group not only had a poor response to sorafenib, but also had a significantly shorter overall survival, suggesting that our RS model can predict treatment response and long-term prognosis at the same time. It is worth noting that DNASE1L3 and ACSL6 of the five characteristic genes have been reported to be associated with the progression and metastasis of HCC (31,32). In particular, DNASE1L3 has been shown to be an independent predictor of prognosis after HCC resection (33). These findings support that our model gene is not only involved in the regulation of sorafenib responsiveness, but also affects the invasive biological behavior of tumors. Therefore, by capturing the inherent biological characteristics of tumor, the model can achieve the dual prediction of treatment response and prognosis.

The following study selected hub modules based on the HCC sensitivity to sorafenib via WGCNA. On the basis of TCGA and GSE109211 datasets, a total of 399 significantly DEGs associated with sorafenib response were obtained. GO enrichment analysis revealed that these genes were associated with cell proliferation and metabolic processes. Study has revealed that sorafenib response is related to a metabolic shift towards aerobic glycolysis (34). Moreover, sorafenib could inhibit the expression of membrane-bound complement regulatory protein, thereby enhancing the antitumor effects of ofatumumab and rituximab in chronic lymphocytic leukaemia (35). Interestingly, KEGG analysis showed that these DEGs were involved in complement and coagulation cascades pathway. In addition, solute carrier (SLC) transporters, which mediate cellular uptake of sorafenib, are critical determinants of drug sensitivity. Our KEGG pathway analysis identified drug metabolism-cytochrome P450 as a significantly enriched pathway, which is closely linked to SLC-mediated drug transport. Although the five signature genes in our model are not direct SLC family members, their association with HCC sensitivity to sorafenib suggests potential indirect regulation of drug uptake. For instance, ACSL6, involved in fatty acid metabolism (36), may affect membrane lipid composition and thus SLC transporter activity. CD68⁺ macrophages, enriched in high-risk groups, could secrete cytokines that downregulate SLC transporters in tumor cells, contributing to resistance (37). Future studies should explore the interaction between the signature genes and SLC family members (e.g., SLC22A1) to clarify their role in sorafenib uptake. Taken together, these results indicated that the metabolic and complement pathways and SLC transporters may be associated with sorafenib response in HCC. However, the roles of these pathways and SLC transporters in sorafenib resistance warranted to be further investigated in HCC in the future.

By comparing the 399 DEGs with module genes screened by WGCNA, 56 sorafenib response-related genes were screened. Finally, five signature genes (DNASE1L3, ACSL6, ACAN, BRSK1 and CD68) were obtained after LASSO analysis. DNASE1L3 has been reported to inhibit HCC progression by weakening glycolysis (38). Expression of DNASE1L3 is an independent prognostic factor of better prognosis in HCC patients after resection (39). ACSL6 belongs to long-chain Acyl-CoA synthetase family that plays a critical role in fatty acid metabolism (40). Upregulation of ACSL6 can promote liver cell apoptosis in the animal model of non-alcoholic fatty liver disease (41). ACSL4 has been suggested to be a predictive biomarker of sorafenib sensitivity in HCC (42), but the role of ACSL6 in HCC has not been elucidated. ACAN is a member of the aggrecan/versican proteoglycan family present in the extracellular matrix of a variety of tissues (43). HCC-derived label-retaining cancer cells downregulate ACAN after treatment with sorafenib (44). BRSK1 is a AMPK-related kinase, which is amplified in liver tumors compared to normal tissues (45). CD68 is one of the most common tumorassociated macrophages markers (46). The infiltration of CD68 positive cells is associated with poor prognosis in patients with HCC (47). Our study for the first time suggested the potential role of these genes in predicting sorafenib response in HCC.

With the five genes, a RS prognostic risk prediction model was constructed. The ROC curves revealed that the 1-, 3- and 5-year AUC values were all above 0.7 in both training and validation sets, suggesting that the model constructed by the five genes presented high predictive precision. Additionally, a nomogram was constructed based on pathologic stage and RS model status, the nomogram also presented good prediction performance in survival prognosis. Overall, we inferred that the five genes can be used as candidate genes to predict sorafenib response in HCC, but the exact role needs to be further elucidated.

According to the RS values, the samples were divided into high- and low-risk groups, and the existed obvious differences on immunity between two groups. Cabrera et al. (48) have reported that in addition to direct action on tumor cells to cause apoptosis and anti-angiogenesis effects, sorafenib can also enhance tumor immunity. The TME, including immune cells, stromal cells, and extracellular matrix, plays a pivotal role in sorafenib response, a factor that is often overlooked in traditional cell culture models. In our study, we explicitly analyzed TME features by evaluating immune cell infiltration and stromal/immune scores, revealing that risk groups differed significantly in immune composition. High-risk groups were enriched with immunosuppressive cells (e.g., Tregs, M0 macrophages) and had lower immune scores, while low-risk groups showed higher infiltration of effector immune cells. These findings align with previous reports that a pro-inflammatory TME enhances sorafenib sensitivity by promoting drug penetration and anti-tumor immunity (49). Conversely, a fibrotic or immunosuppressive TME can limit sorafenib diffusion and blunt its efficacy, which may explain the poor response in high-risk groups (50). Our model integrates TME-related genes (e.g., CD68, a macrophage marker) and immune infiltration data, partially addressing the influence of TME on drug response. Furthermore, Tregs have been reported to exert immunosuppressive effect in HCC patients by suppressing beneficial antitumor immunity (51). Clinically, a high number of Tregs is an indicator of poor prognosis, whereas a low number of Tregs is associated with improved survival (52). In our study, the number of Tregs in the high risk group was significantly more than that in the low risk group, in accordance with the report above. It has been reported that the proportion of immune cells and tumor purity in the TME is associated with the prognosis of HCC, and patients with abundant immune cell infiltration had better response to treatment (53). Therefore, we speculated that the patients in the high risk group may be benefit from therapy. The finding that the high-risk group had a lower sorafenib IC₅₀ value also supported our hypothesis.

In addition, sorafenib exerts dual effects on tumor cells and the immune system through tyrosine kinase inhibition (54). Mechanistically, sorafenib blocks VEGF/VEGFR signaling, which not only inhibits angiogenesis but also reduces VEGF-mediated immunosuppression-VEGF is known to impair dendritic cell maturation and T cell infiltration (55,56). Additionally, sorafenib inhibits RAF/MEK/ERK pathways in tumor cells, thereby downregulating the expression of immunosuppressive molecules such as programmed cell death-ligand 1 (PD-L1) and interleukin-10 (IL-10) (57,58). Our finding that high-risk groups (associated with sorafenib resistance) have enriched Tregs and M2-like macrophages aligns with these mechanisms: ineffective tyrosine kinase inhibition in resistant cells may preserve the immunosuppressive microenvironment. Conversely, effective sorafenib treatment in low-risk groups could relieve immune suppression, promoting anti-tumor immunity.

Although this study focuses on the prediction of sorafenib response, our findings also provide a potential clue for predicting the response to immunotherapy. The close correlation between our risk model and immune cell infiltration suggests its potential utility in predicting response to immunotherapy. The five signature genes are involved in immune regulation, providing a theoretical basis for this inference. CD68, a marker of tumor-associated macrophages (TAMs), is highly expressed in the high-risk group. TAMs promote immunosuppression by secreting cytokines such as IL-10 and transforming growth factor-beta (TGF-β), which is associated with resistance to immune checkpoint inhibitors (ICIs) (59,60). Thus, high CD68 expression in the high-risk group may indicate poor response to ICIs alone, but combination with TAM-depleting agents (e.g., anti-CSF1R) could improve efficacy. Furthermore, Tregs are enriched in the high-risk group, and Tregs-mediated immunosuppression is a major barrier to ICIs efficacy. Studies have shown that patients with high Tregs infiltration benefit more from anti-CTLA-4 therapy, which depletes Tregs (61,62). The ability of our model to identify Tregs-enriched high-risk groups may help select patients suitable for CTLA-4 inhibitors. DNASE1L3, which inhibits glycolysis in HCC (38), may enhance immunotherapy response. Glycolysis by tumor cells consumes glucose, limiting T cell function (63). Low DNASE1L3 expression in the high-risk group may exacerbate this metabolic competition, but combining ICIs with glycolysis inhibitors (e.g., metformin) could reverse it. These findings suggest that our model may have the potential to predict the response to immunotherapy.

However, there are some limitations in this study. First, our proposed model requires validation in HCC cohorts treated with immunotherapy (e.g., atezolizumab plus bevacizumab). Future studies should integrate our model with immunotherapy response data to confirm its predictive value. Second, our study did not directly analyze the role of SLC transporters in sorafenib resistance, which is a known critical factor. The relationship between the five signature genes and SLC-mediated drug uptake remains to be validated experimentally, such as through in vitro assays of sorafenib accumulation in cells with altered expression of these genes. Additionally, additional clinical factors that affecting HCC prognosis, such as alpha-fetoprotein levels, Eastern Cooperative Oncology Group performance status, and Child-Pugh classification, should be further included.

Conclusions

In conclusion, our study established a sorafenib response-related prognostic risk prediction model in HCC based on five signature genes (DNASE1L3, ACSL6, ACAN, BRSK1 and CD68), which had high predictive precision. There were significant correlations between risk groups and immunity. The patients in the high risk group may be benefit from sorafenib treatment. These findings provided predictive biomarkers for sorafenib response in HCC patients, which may contribute to the precision therapy for HCC.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-1020/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-1020/prf

Funding: This study was supported by the Prevention and Treatment of AMR via B7H4-LILRB4 Cascade Regulating Metabolic Reprogramming of B-Cell during Alloantibody Production (No. 82371792).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2024-1020/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. Ferlay J, Soerjomataram I, Dikshit R, et al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer 2015;136:E359-86. [Crossref] [PubMed]
2. Yuan W, Tao R, Huang D, et al. Transcriptomic characterization reveals prognostic molecular signatures of sorafenib resistance in hepatocellular carcinoma. Aging (Albany NY) 2021;13:3969-93. [Crossref] [PubMed]
3. Choi KJ, Baik IH, Ye SK, et al. Molecular Targeted Therapy for Hepatocellular Carcinoma: Present Status and Future Directions. Biol Pharm Bull 2015;38:986-91. [Crossref] [PubMed]
4. Ocana A, Pandiella A, Siu LL, et al. Preclinical development of molecular-targeted agents for cancer. Nat Rev Clin Oncol 2010;8:200-9. [Crossref] [PubMed]
5. Kudo M, Finn RS, Qin S, et al. Lenvatinib versus sorafenib in first-line treatment of patients with unresectable hepatocellular carcinoma: a randomised phase 3 non-inferiority trial. Lancet 2018;391:1163-73. [Crossref] [PubMed]
6. Luo J, Gao B, Lin Z, et al. Efficacy and safety of lenvatinib versus sorafenib in first-line treatment of advanced hepatocellular carcinoma: A meta-analysis. Front Oncol 2022;12:1010726. [Crossref] [PubMed]
7. Forner A, Reig M, Bruix J. Hepatocellular carcinoma. Lancet 2018;391:1301-14. [Crossref] [PubMed]
8. Tang W, Chen Z, Zhang W, et al. The mechanisms of sorafenib resistance in hepatocellular carcinoma: theoretical basis and therapeutic aspects. Signal Transduct Target Ther 2020;5:87. [Crossref] [PubMed]
9. Villanueva A, Llovet JM. Second-line therapies in hepatocellular carcinoma: emergence of resistance to sorafenib. Clin Cancer Res 2012;18:1824-6. [Crossref] [PubMed]
10. Marisi G, Cucchetti A, Ulivi P, et al. Ten years of sorafenib in hepatocellular carcinoma: Are there any predictive and/or prognostic markers? World J Gastroenterol 2018;24:4152-63. [Crossref] [PubMed]
11. Arao T, Ueshima K, Matsumoto K, et al. FGF3/FGF4 amplification and multiple lung metastases in responders to sorafenib in hepatocellular carcinoma. Hepatology 2013;57:1407-15. [Crossref] [PubMed]
12. Horwitz E, Stein I, Andreozzi M, et al. Human and mouse VEGFA-amplified hepatocellular carcinomas are highly sensitive to sorafenib treatment. Cancer Discov 2014;4:730-43. [Crossref] [PubMed]
13. de Visser KE, Joyce JA. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 2023;41:374-403. [Crossref] [PubMed]
14. Hanahan D, Monje M. Cancer hallmarks intersect with neuroscience in the tumor microenvironment. Cancer Cell 2023;41:573-80. [Crossref] [PubMed]
15. Pinyol R, Montal R, Bassaganyas L, et al. Molecular predictors of prevention of recurrence in HCC with sorafenib as adjuvant treatment and prognostic factors in the phase 3 STORM trial. Gut 2019;68:1065-75. [Crossref] [PubMed]
16. Grinchuk OV, Yenamandra SP, Iyer R, 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. [Crossref] [PubMed]
17. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
18. Zhu Y, Meng X, Ruan X, et al. Characterization of Neoantigen Load Subgroups in Gynecologic and Breast Cancers. Front Bioeng Biotechnol 2020;8:702. [Crossref] [PubMed]
19. Wang P, Wang Y, Hang B, et al. A novel gene expression-based prognostic scoring system to predict survival in gastric cancer. Oncotarget 2016;7:55343-51. [Crossref] [PubMed]
20. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
21. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
22. Goeman JJ. L1 penalized estimation in the Cox proportional hazards model. Biom J 2010;52:70-84. [Crossref] [PubMed]
23. Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med 1997;16:385-95. [Crossref] [PubMed]
24. Harrell FE Jr, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med 1996;15:361-87. [Crossref] [PubMed]
25. Shan S, Chen W, Jia JD. Transcriptome Analysis Revealed a Highly Connected Gene Module Associated With Cirrhosis to Hepatocellular Carcinoma Development. Front Genet 2019;10:305. [Crossref] [PubMed]
26. Chen B, Khodadoust MS, Liu CL, et al. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol 2018;1711:243-59. [Crossref] [PubMed]
27. Hu D, Zhou M, Zhu X. Deciphering Immune-Associated Genes to Predict Survival in Clear Cell Renal Cell Cancer. Biomed Res Int 2019;2019:2506843. [Crossref] [PubMed]
28. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
29. Heimbach JK, Kulik LM, Finn RS, et al. AASLD guidelines for the treatment of hepatocellular carcinoma. Hepatology 2018;67:358-80. [Crossref] [PubMed]
30. Brown ZJ, Tsilimigras DI, Ruff SM, et al. Management of Hepatocellular Carcinoma: A Review. JAMA Surg 2023;158:410-20. [Crossref] [PubMed]
31. Zhang H, Liu R, Sun L, et al. The Effect of Ferroptosis-Related Genes on Prognosis and Tumor Mutational Burden in Hepatocellular Carcinoma. J Oncol 2021;2021:7391560. [Crossref] [PubMed]
32. Wang J, Chen Y, Xu Y, et al. DNASE1L3-mediated PANoptosis enhances the efficacy of combination therapy for advanced hepatocellular carcinoma. Theranostics 2024;14:6798-817. [Crossref] [PubMed]
33. Deng Z, Xiao M, Du D, et al. DNASE1L3 as a Prognostic Biomarker Associated with Immune Cell Infiltration in Cancer. Onco Targets Ther 2021;14:2003-17. [Crossref] [PubMed]
34. Castello A, Rimassa L, Personeni N, et al. Metabolic Switch in Hepatocellular Carcinoma Patients Treated with Sorafenib: a Proof-of-Concept Trial. Mol Imaging Biol 2020;22:1446-54. [Crossref] [PubMed]
35. Dwojak M, Bobrowicz M, Bil J, et al. Sorafenib improves rituximab and ofatumumab efficacy by decreasing the expression of complement regulatory proteins. Blood Cancer J 2015;5:e300. [Crossref] [PubMed]
36. Li Z, Xiang YJ, Zou ZC, et al. Multi-omics analysis and experimental verification reveal testicular fatty acid metabolism disorder in non-obstructive azoospermia. Zool Res 2025;46:177-92. [Crossref] [PubMed]
37. Ansari RE, Craze ML, Althobiti M, et al. Enhanced glutamine uptake influences composition of immune cell infiltrates in breast cancer. Br J Cancer 2020;122:94-101. [Crossref] [PubMed]
38. Xiao Y, Yang K, Liu P, et al. Deoxyribonuclease 1-like 3 Inhibits Hepatocellular Carcinoma Progression by Inducing Apoptosis and Reprogramming Glucose Metabolism. Int J Biol Sci 2022;18:82-95. [Crossref] [PubMed]
39. Wang S, Ma H, Li X, et al. DNASE1L3 as an indicator of favorable survival in hepatocellular carcinoma patients following resection. Aging (Albany NY) 2020;12:1171-85. [Crossref] [PubMed]
40. Yan S, Yang XF, Liu HL, et al. Long-chain acyl-CoA synthetase in fatty acid metabolism involved in liver and other diseases: an update. World J Gastroenterol 2015;21:3492-8. [Crossref] [PubMed]
41. Xu C, Wang G, Hao Y, et al. Correlation analysis between gene expression profile of rat liver tissues and high-fat emulsion-induced nonalcoholic fatty liver. Dig Dis Sci 2011;56:2299-308. [Crossref] [PubMed]
42. Feng J, Lu PZ, Zhu GZ, et al. ACSL4 is a predictive biomarker of sorafenib sensitivity in hepatocellular carcinoma. Acta Pharmacol Sin 2021;42:160-70. [Crossref] [PubMed]
43. Vafaeie F, Nomiri S, Ranjbaran J, et al. ACAN, MDFI, and CHST1 as Candidate Genes in Gastric Cancer: A Comprehensive Insilco Analysis. Asian Pac J Cancer Prev 2022;23:683-94. [Crossref] [PubMed]
44. Xin HW, Ambe CM, Hari DM, et al. Label-retaining liver cancer cells are relatively resistant to sorafenib. Gut 2013;62:1777-86. [Crossref] [PubMed]
45. Nguyen K, Hebert K, McConnell E, et al. LKB1 signaling and patient survival outcomes in hepatocellular carcinoma. Pharmacol Res 2023;192:106757. [Crossref] [PubMed]
46. Pulford KA, Rigney EM, Micklem KJ, et al. KP1: a new monoclonal antibody that detects a monocyte/macrophage associated antigen in routinely processed tissue sections. J Clin Pathol 1989;42:414-21. [Crossref] [PubMed]
47. Kong LQ, Zhu XD, Xu HX, et al. The clinical significance of the CD163+ and CD68+ macrophages in patients with hepatocellular carcinoma. PLoS One 2013;8:e59771. [Crossref] [PubMed]
48. Cabrera R, Ararat M, Xu Y, et al. Immune modulation of effector CD4+ and regulatory T cell function by sorafenib in patients with hepatocellular carcinoma. Cancer Immunol Immunother 2013;62:737-46. [Crossref] [PubMed]
49. Zhu J, Han T, Zhao S, et al. Computational Characterizing Necroptosis Reveals Implications for Immune Infiltration and Immunotherapy of Hepatocellular Carcinoma. Front Oncol 2022;12:933210. [Crossref] [PubMed]
50. Li C, Xiong L, Yang Y, et al. Sorafenib enhanced the function of myeloid-derived suppressor cells in hepatocellular carcinoma by facilitating PPARα-mediated fatty acid oxidation. Mol Cancer 2025;24:34. [Crossref] [PubMed]
51. Zahran AM, Nafady-Hego H, Mansor SG, et al. Increased frequency and FOXP3 expression of human CD8(+)CD25(High+) T lymphocytes and its relation to CD4 regulatory T cells in patients with hepatocellular carcinoma. Hum Immunol 2019;80:510-6. [Crossref] [PubMed]
52. Zhou J, Ding T, Pan W, et al. Increased intratumoral regulatory T cells are related to intratumoral macrophages and poor prognosis in hepatocellular carcinoma patients. Int J Cancer 2009;125:1640-8. [Crossref] [PubMed]
53. Ge PL, Li SF, Wang WW, et al. Prognostic values of immune scores and immune microenvironment-related genes for hepatocellular carcinoma. Aging (Albany NY) 2020;12:5479-99. [Crossref] [PubMed]
54. Eun JW, Yoon JH, Ahn HR, et al. Cancer-associated fibroblast-derived secreted phosphoprotein 1 contributes to resistance of hepatocellular carcinoma to sorafenib and lenvatinib. Cancer Commun (Lond) 2023;43:455-79. [Crossref] [PubMed]
55. Kalathil SG, Hutson A, Barbi J, et al. Augmentation of IFN-γ+ CD8+ T cell responses correlates with survival of HCC patients on sorafenib therapy. JCI Insight 2019;4:e130116. [Crossref] [PubMed]
56. Hung SW, Zhang R, Tan Z, et al. Pharmaceuticals targeting signaling pathways of endometriosis as potential new medical treatment: A review. Med Res Rev 2021;41:2489-564. [Crossref] [PubMed]
57. Liang Y, Chen J, Yu Q, et al. Phosphorylated ERK is a potential prognostic biomarker for Sorafenib response in hepatocellular carcinoma. Cancer Med 2017;6:2787-95. [Crossref] [PubMed]
58. Tang F, Li S, Liu D, et al. Sorafenib sensitizes melanoma cells to vemurafenib through ferroptosis. Transl Cancer Res 2020;9:1584-93. [Crossref] [PubMed]
59. Pan Y, Yu Y, Wang X, et al. Tumor-Associated Macrophages in Tumor Immunity. Front Immunol 2020;11:583084. [Crossref] [PubMed]
60. Xiang X, Wang J, Lu D, et al. Targeting tumor-associated macrophages to synergize tumor immunotherapy. Signal Transduct Target Ther 2021;6:75. [Crossref] [PubMed]
61. Régnier P, Vetillard M, Bansard A, et al. FLT3L-dependent dendritic cells control tumor immunity by modulating Treg and NK cell homeostasis. Cell Rep Med 2023;4:101256. [Crossref] [PubMed]
62. Guo C, Dai X, Du Y, et al. Preclinical development of a novel CCR8/CTLA-4 bispecific antibody for cancer treatment by disrupting CTLA-4 signaling on CD8 T cells and specifically depleting tumor-resident Tregs. Cancer Immunol Immunother 2024;73:210. [Crossref] [PubMed]
63. Geiger R, Rieckmann JC, Wolf T, et al. L-Arginine Modulates T Cell Metabolism and Enhances Survival and Anti-tumor Activity. Cell 2016;167:829-842.e13. [Crossref] [PubMed]