Introduction

Hepatocellular carcinoma (HCC) ranks as the sixth most common malignancy and the fourth leading cause of cancer-related death worldwide.1 Approximately 80% of cases are diagnosed at intermediate or advanced stages.2 Transarterial chemoembolization (TACE) is recognized as the primary treatment modality for intermediate-stage HCC, where it targetedly delivers chemotherapeutics while blocking the tumor’s blood supply.3 In addition, TACE can be applied in early-stage HCC who are waiting for a liver transplant, those not feasible for curative treatment, patients receiving adjuvant therapy for post-curative treatment, recurrent lesions after curative surgical resection or ablation, and advanced HCC combining with systemic therapy.4–6 Nevertheless, the therapeutic efficacy of TACE varies among individuals, leading to the emergence of the concept of TACE refractoriness/failure.7 The exact underlying mechanism of TACE failure remains elusive, with several potential explanations, mainly including VEGF upregulation-induced angiogenesis and neovascularization and resistance to chemotherapeutic drugs after consecutive TACE therapies.8 There is a pressing need to identify biomarkers that can predict TACE response preoperatively and help select patients who are most likely to benefit from TACE.

Due to rapid advancements in RNA sequencing, an increasing number of signatures based on transcriptional information have been developed to comprehensively and deeply elucidate the molecular mechanism underlying tumor onset or progression, establish correlations with tumor immune microenvironment characteristics and offer prognostic insights.9,10 A series of TACE refractoriness prediction scores have been established by integrating clinical variables, such as tumor count, “up-to-seven” criteria application status, and bilobular invasion status.11,12 Furthermore, a gene signature has been constructed to predict TACE failure and the prognosis of HCC for patients undergoing TACE therapy.13,14 However, the widespread use of high-throughput sequencing in clinical practice is hindered by its relatively high cost and invasive procedures.

Contrast-enhanced computed tomography (CECT) or magnetic resonance imaging (MRI) are standardized imaging modalities for HCC diagnosis, prognostication and follow-up.15 Besides visible or measurable parameters in CT or MRI, visible imaging data can be converted into high-throughput and precise features through radiomics approach, and it has been extensively used in various cancer types for diagnosis and clinical outcome prediction.16,17 Nevertheless, the biological significance underlying the radiomics features remains to be further clarified. In recent years, an emerging method known as radiogenomics can link imaging phenotypes to tumor biological characteristics, thereby enabling noninvasive “digital biopsy”.18–21

Therefore, this study aimed to establish a transcriptomic biomarker to predict the response to TACE, which helps select best HCC candidates for TACE treatment. In addition, the correlations between gene signature with CT-based radiomics features and tumor immune microenvironment characteristics in patients with HCC were explored.

Materials and Methods Data Collection

Gene expression data from the GSE104580 dataset were downloaded from the Gene Expression Omnibus (GEO) database, containing 81 tissue RNA samples from TACE_sensitive patients and 66 tissue RNA samples from TACE_resistant patients. All samples were derived from tumor biopsies of patients with HCC prior to TACE. The response to TACE within 3 months (after the first or second TACE session) was assessed by extramural reviewers via the modified Response Evaluation Criteria in Solid Tumors (mRECIST). In this study, the tumor response was the primary endpoint. Patients with a complete response or a partial response were defined as having an objective response to TACE, whereas patients with stable disease or progressive disease were determined to be TACE-resistant.22

The transcriptomic, mutational and clinicopathological data of another HCC cohort were obtained from The Cancer Genome Atlas Liver Hepatocellular Carcinoma dataset (containing 374 hCC samples and 50 control samples) (TCGA-LIHC; https://portal.gdc.cancer.gov/). The International Cancer Genome Consortium-Liver Cancer-Riken-Japan (ICGC-LIRI-JP) dataset was used to assess the prognostic value of the gene signature. The expression profiles were normalized and background adjusted using the “limma” R package. The corresponding computed tomography (CT) images of the TCGA-LIHC cohort were retrieved from The Cancer Imaging Archive (TCIA; https://www.cancerimagingarchive.net/). This study is reported in accordance with the Standards for the Reporting of Diagnostic Accuracy Studies (STARD) criteria.23

This study was approved by the Ethics Committee of our center and adhered to the ethical principles of the Declaration of Helsinki. The informed consent was waived due to the nature of retrospective study.

Hub Gene Selection and Signature Construction

First, differential gene expression analysis between TACE-sensitive and TACE-resistant samples was conducted with the “limma” package, with false discovery rate (FDR) values < 0.05 and |logFC| >1.0. A heatmap and volcano plot were generated to visualize the differentially expressed genes (DEGs).

Subsequently, the machine learning algorithms least absolute shrinkage and selection operator (LAASO), random forest (RF) and support vector machine-recursive feature elimination (SVM-RFE) were utilized to screen hub genes from the DEGs. A nomogram was constructed based on the overlapping results from three machine learning methods to visualize the multivariate logistic regression analysis. The gene signature was generated with the formula: TACE Failure Score=Constant+ (Expression of gene 1*Coefficient of gene 1) +(Expression of gene 2*Coefficient of gene 2) +……+ (Expression of gene i*Coefficient of gene i). The prediction efficiency was evaluated using receiver operating characteristic (ROC) analysis, calibration, and decision curve analysis (DCA).

Gene Co-Expression Network Analysis

The “WGCNA” package was applied to establish the coexpression network between TACE-Nonresponse and TACEResponse samples. The integrity of data was checked with the “goodSampleGenes”function. Then, pearson’s correlation analysis of all pairs of genes was used to establish an adjacency matrix, which was used to construct a scale-free co-expression network. The adjacency matrix was transformed to a topological overlap matrix (TOM) that could quantify the network connectivity of a gene. Following the calculation of module eigengene (ME) and merging similar modules in the clustering tree based on ME, a hierarchical clustering dendrogram was generated. Gene significance (GS), which was the mediator p-value (GS = lgP) for each gene, represented the degree of linear correlation between gene expression and clinical states. Module Membership (MM) > 0.5 and GS > 0.3 were set as the threshold to screen hub genes in each module.

Imaging Processing and Analysis

The inclusion criteria for CT imaging (CTI) were as follows: (1) preoperative imaging; (2) contrast-enhanced scans; and (3) matched clinical and genomic information. Overall, thirty-six patients were included in the subsequent analysis.

Lesions were identified and manually delineated via ITK SNAP software (open source, www.itk-snap.org) on each slice of CECT images by a radiologist with 10 years of abdominal radiology experience, and the observations were validated by a senior radiologist with 15 years of experience.

Radiomic features were extracted with the “PyRadiomics” package (https://www.radiomics.io/pyradiomics.html) and included three sets: (1) first-order statistics; (2) shape and size; and (3) high-order statistical textural information.

The interobserver reproducibility of radiomic feature extraction was assessed using intraclass correlation coefficients (ICCs), with 0.81–1.00 indicating almost perfect agreement, 0.61–0.80 indicating substantial agreement, and 0.41–0.60 indicating moderate agreement. Features with an ICC > 0.75 were retained for further analysis. Subsequently, all features were normalized using the min–max approach to achieve a unit norm. T tests and LASSO analyses were used to identify features strongly correlated with the TACE response. Next, the selected radiomic features were integrated into a Rad-score model using a logistic regression formula as follows: Rad score=Constant+ (feature 1*Coefficient of feature 1) + (feature 2*Coefficient of feature 2) + …… + (feature i*Coefficient of feature i). ROC analysis was conducted to evaluate the diagnostic value for TACE failure. In addition, the correlations between the Rad-score and hub gene expression and between the Rad-score and TFS score were examined by Spearman analysis.

Patient Enrollment and Follow-Up

For external validation, patients who were initially treated with TACE at our center from January 2019 to June 2019 were retrospectively reviewed. The inclusion criteria were as follows: 1) age≥18 years, 2) histologically or clinically diagnosed with HCC, 3) BCLC stage B or C, 3) preserved liver function (Child–Pugh A grade), 4) Eastern Cooperative Oncology Performance Status (ECOG-PS) 0–1, and 4) contrast-enhanced CT examination performed within 2 weeks before TACE treatment. Patients who 1) had received other treatments, such as transplantation, surgery, ablation, or systemic therapies before the initial TACE treatment; 2) had a history of other cancers; 3) lacked necessary clinical information on demographic characteristics, laboratory examinations and tumor characteristics; or 4) lacked adequate image quality for further lesion delineation were excluded.

Contrast-enhanced CT or MRI was performed within 4–6 weeks after the operation for follow-up and response evaluation based on the modified Response Evaluation Criteria in Solid Tumors (mRECIST) criteria. Responses to the first session of TACE were divided into complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). Patients with a CR or PR were defined as responders; those with PD and SD were defined as nonresponders. In addition, the criteria for defining TACE-sensitive and TACE-resistant patients in our cohort was consistent with the GSE104580 dataset.

Associations of the Gene Signature with Clinical Characteristics, Tumor Mutation Burden and Survival

Differences in clinical characteristics, including age, sex, histologic grade, TNM stage, and immune subtypes, between risk groups were measured. To determine the prognostic value of the gene signature, Kaplan‒Meier analysis was performed with the “survival” and “survminer” R packages. Somatic mutation profiles were downloaded from the General Genomic Data (GDC) website. The “maftools” package was used to evaluate tumor mutation burden (TMB).

Enrichment Analysis

The package “clusterProfiler” was utilized to perform Gene Ontology (GO) enrichment analysis to identify top five types of significantly enriched GO terms (p<0.05) for clarifying the functional roles of DEGs. Gene set enrichment analysis (GSEA) was performed to identify underlying enriched pathways among the high- and low-risk groups through Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (https://www.gsea-msigdb.org/gsea/index.jsp). P values <0.05, |NES| >1 and FDR q < 0.25 were considered to indicate statistical significance.

Expressions of Key Genes in HCC Cell Lines and Single-Cell Level

Expression levels of key genes in 24 kinds of HCC cell lines were verified utilizing the Human Protein Atlas (HPA) online database. Additionally, the expression of key genes at single-cell level was further investigated through analyzing the data from HCCDB v2.0 repository (available at https://lifeome.net/database/hccdb2).

Immune Cell Infiltration Analysis and Immunotherapy Response Evaluation

The extent of immune cell infiltration in all samples was assessed using the CIBERSORT algorithm and single-sample gene set enrichment analysis (ssGSEA). CIBERSORT is a computational algorithm for transforming a gene expression matrix to reveal the levels of 22 types of infiltrating immune cells, with the threshold for statistical significance set at a p value <0.01. The ssGSEA algorithm was based on 29 immune gene sets, including genes related to different immune cell types, functions, pathways, and checkpoints. The R packages “GSVA”, “GSEABase”, and “limma” were used to evaluate the immunologic characteristics of the samples. ESTIMATE analysis was used to compute the immune, stromal and ESTIMATE scores.

The immunotherapy efficacy in patients in the TCGA-LIHC cohort was evaluated according to the immunophenoscore (IPS) obtained from the Cancer Immunome Atlas (http://tcia.at/). Additionally, the IMvigor210 cohort was used to evaluate the ability of the signature to predict immunotherapy response. Comparisons of immune checkpoint molecule expression between the two risk groups were also conducted. The association between the gene signature and the response to sorafenib treatment was analyzed in the GSE109211 dataset, which included data from 21 patients who were sorafenib responders and 46 who were nonresponders.

Durg Sensitivity Analysis

The “oncoPredict” R package was used to detect differences in sensitivity to commonly used chemotherapy and molecular targeted drugs between the two risk groups.

Results TACE Failure Gene Signature Establishment

The study flowchart is presented in Figure 1. A total of 276 genes were identified by differential expression analysis (Figure 2a and b). A Venn diagram was generated to illustrate the five hub genes selected by the three machine learning algorithms (Figure 2c–e), namely, ADH1C, CXCL11, EMCN, SPARCL1, and LIN28B (Figure 2f). Genes selected by machine learning algorithms were listed in Supplementary Table S1. Following WGCNA analysis, a total of 16 modules were determined with the dynamic tree cutting method, and each color represented one module. Nine hub modules were identified through module-trait correlation analysis, among which the turquoise module was considered most significant. The above 5 key genes were present in significant modules (Figure 2g–l).

Figure 1 The flowchart of this study.

Figure 2 Identification of hub genes. (a) Heatmap of differentially expressed genes, (b) Volcano plot for visualization of upregulated and downregulated genes between TACE responders and TACE nonresponders, (c) Support vector machine recursive feature elimination (SVM-RFE) results, (d) Least absolute shrinkage and selection operator (LASSO) regression analysis results, (e) Random forest (RF) analysis results, (f) Venn diagram for intersection of the three machine learning methods, (g) Analysis of network topology for various soft thresholding powers. (h) clustering dendrogram of mRNAs, with dissimilarity based on topological overlap, together with assigned module colors. (i) Heatmap showing the correlation between modules and clinical phenotype. (j–l) Scatter plots of module eigengenes in turquoise module, yellow and brown module.

A logistic regression algorithm was used to fit the prediction model and generate the TACE failure signature (TFS). A nomogram was constructed to visualize the results of the multivariate logistic regression analysis (Figure 3a). The expression levels of five hub genes and the TFS score were compared between TACE responders and TACE nonresponders (Figure 3b and c). The area under the ROC curve for the gene signature was 0.916 (Figure 3d). Both the calibration curve and decision curve demonstrated satisfactory prediction accuracy and net clinical benefits (Figure 3e and f). The TFS formula was as follows:

Figure 3 TACE failure signature (TFS) construction in Gene Expression Omnibus (GEO) cohort. (a) Nomogram for predicting TACE response based on the coefficients of multivariate logistic regression analysis, (b) Gene expression levels of five hub genes, (c) Comparison of TFS scores between TACE responders and TACE nonresponders, (d) Receiver operating characteristic (ROC) analysis showing model discrimination ability, (e) Calibration curve of the prediction model, (f) Decision curve of the model. *p<0.05, **p<0.01, ***p<0.001.

TFS score=14.7795+ADH1C expression*(−0.3420) + SPARCL1 expression*(−0.9699) +EMCN expression*(−0.2586) +CXCL11 expression*(−0.1831) +LIN28B expression*(0.3229).

Expressions of Key Genes

Expressions of five hub genes in different HCC cell lines exhibits heterogeneity, which may be attributed to multiple factors, including the inherent heterogeneity of liver cancer, genomic and epigenetic variations between cell lines, differences in signaling pathway activity, and the state of cellular differentiation (Supplementary Figure S1). As depicted in UMAP, and t-SNE plots derived from HCCDB single-cell data, we identified the expression patterns of five key genes across nine major cell clusters. These clusters include B cells, cholangiocytes, endothelial cells, hepatocytes, malignant cells, myeloid cells, NK/T cells, plasma B cells, and stromal cells (Supplementary Figure S2). The results revealed distinct expression profiles for each gene within these cell types, highlighting potential roles in tumor microenvironment dynamics and cellular interactions. This detailed expression mapping underscores the heterogeneous landscape of HCC at the single-cell level, providing insights into the functional relevance of these genes in various cellular contexts.

Enrichment Analysis and Correlation Between the TFS Score and Clinical Characteristics

GO analysis was conducted for up-regulated genes between high-risk and low-risk groups to determine the functional differences. Similar to GO analysis of up-regulated genes between TACE-resistant and TACE-sensitive samples in primary GSE104580 cohort, pathways of cellular response to xenobiotic stimulus, response to xenobiotic stimulus, and xenobiotic metabolic process were functionally enriched (Figure 4a and b).

Figure 4 Enrichment analysis and validation of signature in The Cancer Genome Atlas (TCGA) cohort. (a and b) Gene ontology (GO) analysis and (c and d) Gene set enrichment analysis (GSEA) between the high-risk and low-risk groups. (e and f) Kaplan‒Meier curves of TFS in the TCGA-LIHC cohort and ICGC-LIJI-JP cohort, (g) Distribution patterns of clinical characteristics in different risk groups, (h) Comparisons of tumor mutation burden (TMB) between the two risk groups and (i) Correlation results between the TFS score and TMB.

GSEA revealed that the genes upregulated in the high-risk group were enriched in pathways related to the cell cycle, glycosphingolipid biosynthesis (lacto-and-neolacto-series), neuroactive ligand‒receptor interaction, ribosome, and sphingolipid metabolism (Figure 3m). Genes upregulated in the low-risk group were enriched in pathways related to complement and coagulation cascades, drug metabolism (cytochrome P450), fatty acid metabolism, retinol metabolism, and steroid hormone biosynthesis (Figure 4c and d).

Survival in the low-risk group (stratified by the median TFS) was superior to that in the high-risk group in both the TCGA-LIHC cohort and the ICGC-LIRI-JP cohort (Figure 4e and f). Six immune subtypes were identified based on the expression profiles across 33 diverse cancer types from TCGA,24 including the wound healing (C1), IFN-γ-dominant (C2), inflammatory (C3), lymphocyte-depleted (C4), immunologically quiet (C5), and TGF beta-dominant (C6) subtypes. The chi-square test results indicated a significant difference in the distribution of immune subtypes between the two groups (Figure 4g). There was no significant difference in the TMB between the two risk groups (Figure 4h), and there was no correlation between the TFS score and TMB (Figure 4i).

Radiogenomic Analysis

A total of 1502 radiomics features were extracted from the CECT images (Figure 5a). According to the ICC analysis and univariate analysis, these features were reduced to twenty-seven. Six features were ultimately selected through LASSO regression analysis (Figure 5b). The weight of each feature is displayed in Figure 5d.

Figure 5 The flowchart and results of radiogenomics analysis. (a) Imaging data acquisition from The Cancer Imaging Archive (TCIA) dataset, tumor segmentation and radiomics feature extraction, (b) Feature selection process of LASSO regression analysis, (c) Correlations between six radiomics features, (d) Weight plot of each feature, (e) Rad-score difference between different risk groups, (f) Receiver operating characteristic (ROC) curve of the Rad-score in predicting the risk of TACE failure, (g) Correlations between gene expression, risk score and radiomics feature value and the Rad-score, (h) Enrolled patients for validation of radiomics model performance, (i) A representative case of unresectable hepatocellular carcinoma (HCC) treated with TACE with partial response, and (j) ROC curve of the combined model (Rad-score+ECOG PS+tumor size) for evaluating the response to TACE in the external validation cohort.

The radiomics signature was established with the following formula:

Rad-score= 0.66+ (original_shape_Elongation)*0.37+(log-sigma-1-0-mm-3D_glszm_GrayLevelNonUniformityNormalized)* 0.19+(wavelet-LLH_glszm_GrayLevelNonUniformityNormalized)*(−0.39)+(wavelet-LHL_glszm_GrayLevelNonUniformityNormalized)*(−0.16)+ wavelet-HLH_glszm_SmallAreaLowGrayLevelEmphasis*(−0.13)+ wavelet-HHL-gldm-SmallDependenceEmphasis*(0.11)+ (wavelet-HHH_ngtdm_Contrast)*(−0.26).

As shown in the boxplot, the Rad-score was greater in the low-risk group than in the high-risk group (Figure 5e). The AUC of the Rad-score for differentiating patients at high risk from those at low risk of TACE failure was 0.844 (Figure 5f). Furthermore, significant correlations were observed between the TFS score and the Rad-score (Figure 5g).

External Validation

Finally, a total of 39 patients were enrolled in the validation cohort (Figure 5h). Comparisons of demographic and clinical characteristics between the two groups in the validation cohort are listed in Table 1. The TACE responders had a greater percentage of patients with an ECOG PS of 0 than did the TACE nonresponders (73.7% vs 40.0%, p=0.034). In addition, the percentage of patients with a maximal tumor diameter ≥10 cm was greater in the TACE-nonresponders group (35.0% vs 10.5%, p=0.075). There were no significant differences in the other parameters between the two groups (Table 1).

Table 1 Demographic and Clinical Characteristics of Patients in Validation Cohort

The arterial phases of the pretreatment CT images were collected. Following the same imaging protocol used for the TCIA-LIHC dataset, the tumor lesions were delineated, and radiomics features were then extracted from the CT images of our validation cohort. The values of six features screened in the TCIA-LIHC cohort were used to calculate the Rad-scores of the validation cohort based on the established Rad-score formula. Then, a ROC curve was plotted to evaluate the discrimination ability of the Rad-score in the validation cohort. The AUC of the Rad-score for predicting the response to initial TACE in our cohort was 0.676 (95% CI, 0.503–0.850), confirming the predictive value of the Rad-score. Then, the ECOG-PS score, tumor size and the Rad-score were subjected to multivariable logistic regression analysis. After that, a combined model was developed and displayed satisfactory performance in predicting TACE treatment response (AUC=0.756, 95% CI, 0.589–0.880) (Figure 5j).

Relationship Between TFS and Immune Microenvironment Profiles

CIBERSORT analysis demonstrated that the proportion of T_cells_CD4_memory_resting was greater in the low-risk group, while the high-risk group had a greater fraction of M0 macrophages (Figure 6a). Additionally, we explored the correlations between the TFS score and immune cell counts. The results showed that the TFS score was positively correlated with the abundance of M0 macrophages and negatively correlated with the abundance of CD4 memory resting T cells (Figure 6b). The enrichment scores of 29 kinds of immune cells, along with their corresponding immune functions and pathways, were quantified with the ssGSEA method. In comparison to the high-risk group, the low-risk group had greater proportions of B cells, mast cells, and NK cells. Various antigen presentation processes, including APC coinhibition, cytolytic activity, and the type 2 IFN response, were found to be weaker in the high-risk group (Figure 6c). The high-risk group had lower stromal scores and ESTIMATE scores than did the low-risk group (Figure 6d). All of the above results suggested an immunosuppressive microenvironment in the high-risk group.

Figure 6 Tumor immune microenvironment characteristics and response to immunotherapy and targeted therapy between different risk groups. (a) Calculation of the number of infiltrating immune cells between the two groups; (b) Correlations between hub gene expression, gene signature score and immune cells; (c) Single sample Gene Set Enrichment Analysis (ssGSEA) score and (d) Stromal score, immune score, and ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) score comparisons between the high- and low-risk groups; (e) Correlations between the TACE failure signature (TFS) and immune checkpoint molecules; (f) Comparisons of immune checkpoint expression levels between the high-risk and low-risk groups; Immune cell Proportion Score (IPS) of patients treated with Programmed Death-1 (PD-1) (g), Programmed death ligand 1(PD-L1) (h), and cytotoxic T Lymphocyte-associated Antigen-4 (CTLA-4) (i) Immune checkpoint inhibitors; (j) Distribution of the Complete Response (CR)/Partial Response (PR) between the two groups, comparisons of the (k) TFS risk score and (l) survival between the risk groups; (m) Comparison of the TFS score between the sorafenib responders and sorafenib nonresponders in the GSE109211 cohort; (n) Receiver operating characteristic (ROC) curve of the TFS score for predicting the sorafenib response. *p<0.05, **p<0.01, ***p<0.001.

Immune Checkpoint Expression and the Response to Immunotherapy and Targeted Drugs

TFS signature genes and the TFS score were significantly correlated with immune checkpoint gene expression (Figure 6e). The high-risk group exhibited increased expression levels of TNFRSF4, CD276, LAG3, LGALS9, TNFRSF14, TNFRSF18, CTLA-4, PDCD1 (PD-1), and TNFSF15, while the low-risk group exhibited increased expression of NRP1, IDO1, CD28, CD40LG, IDO2, PDCD1LG2, CD44, and CD200 (Figure 6f).

As indicated by the IPS, the low-risk group showed better immunotherapeutic responses to PD-1 and CTLA-4 inhibitors (Figure 6g–i). In the IMvigor210 cohort, the proportion of patients who achieved a CR/PR response in the high-risk group was slightly but not significantly greater than that in the low-risk group (Figure 6j). In addition, the TFS of patients in the CR/PR group was greater than that of patients in the SD/PD group, which may indicate that patients at high risk of TACE failure benefit more from immunotherapy (Figure 6k). KM analysis revealed that high-risk patients receiving immunotherapy had shorter survival than low-risk patients (Figure 6l).

Regarding the sorafenib response, significantly greater TFS was found in sorafenib responders than in sorafenib nonresponders (Figure 6m). The AUC for predicting sorafenib response was 0.899 (Figure 6n).

Drug Sensitivity Analysis

Compared to the low-risk group, the high-risk group was less sensitive to 5−fluorouracil, BPD−00008900, dasatinib, fulvestrant, GDC0810, GSK2606414, I−BRD9, lapatinib, MIRA−1, pevonedistat, pictilisib, SCH772984, taselisib, temozolomide and VX−11e and more sensitive to axitinib, AZ960, AZD1208, AZD1332, AZD2014, AZD6482, AZD8055, BMS−345541, BMS−754807, dihydrorotenone, doramapimod, GSK269962A, IGF1R_3801, IRAK4_4710, JAK_8517, JAK1_8709, JQ1, KRAS (G12C) inhibitor−12, LGK974, mitoxantrone, niraparib, obatoclax mesylate, OF−1, PAK_5339, picolinic acid, ribociclib, ruxolitinib, SB216763, sinularin, TAF1_5496, topotecan and ZM447439 (Figure 7a–u).

Figure 7 Chemotherapeutic drug sensitivity analyses between the high and low-risk groups via the “oncoPredict” package (a-u).

Discussion

In this study, a model was constructed and found to act as a biomarker for estimating TACE failure risk according to transcriptional sequencing data. Ultimately, five genes, namely, ADH1C, CXCL11, EMCN, SPARCL1, and LIN28B, were integrated into the TACE failure signature (TFS). Additionally, it showed connections with radiomics features on pretreatment imaging, tumor immune microenvironment characteristics, and the efficacy of immune-targeted therapy in HCC.

TFS-Related Genes for Evaluating TACE Efficacy

Previous studies have provided insights into how certain genes influence tumor biological behaviour and treatment responses in HCC. For instance, ADH1C has been shown to be significantly downregulated in HCC cells, with experimental data revealing that its silencing promotes the proliferation and migration of HCC cell lines.25 This aligns with the observation that low ADH1C expression is associated with pathways driving tumorigenesis, while higher ADH1C levels correlate with a favorable prognosis in HCC.26,27 The finding reinforces the results of our study, where patients with higher ADH1C expression exhibited a better response to TACE, suggesting that ADH1C may serve as a potential biomarker for TACE responsiveness. An in vivo study by Lau et al revealed that SPARCL1 collaborates with SPARC to mitigate angiogenesis and strongly inhibits tumor growth in HCC cells.28 Similar to the current study, SPARCL1 was included in another gene signature predictive of TACE response (AQP1, FABP4, HERC6, LOX, PEG10, S100A8, SPARCL1, TIAM1, TSPAN8, and TYRO3) that was developed using a machine learning method.13 C-X-C motif chemokine ligand 11 (CXCL11) is a protein-coding gene that is associated with ethmoid sinus adenocarcinoma and ethmoid sinus cancer. Among its related pathways are MIF-mediated glucocorticoid regulation and GPCR downstream signaling. CXCL11 in cancer-associated fibroblasts (CAFs) from HCC was shown to promote migration and tumor metastasis through the circUBAP2/miR-4756/IFIT1/3 and LINC00152/miR-205-5p/CXCL11 axes.29,30 It has been reported that the CXCL9/CXCL10/CXCL11/CXCR3 axis regulates immune cell migration, differentiation, and activation, potentially affecting the response to immunotherapy for cancers.31 Therefore, CXCL11’s association with metastatic potential and immune modulation may indicate its importance in determining the success of TACE treatments. In addition, EMCN, a mucin-like sialoglycoprotein, interferes with the assembly of focal adhesion complexes and inhibits cell-extracellular matrix interactions, which are critical processes in cancer progression.32 While EMCN is overexpressed in HCC compared to adjacent normal liver tissue samples in GSE25097 and GSE14520 datasets,33 qRT-PCR assays have reported that the expression of EMCN was low in HCC, and this phenotype indicated a favorable prognosis, underscoring its potential role as a prognostic marker for TACE procedures.34 In our model development cohort, only LIN28B was upregulated in TACE nonresponders. LIN28B upregulation could counteract the anti-proliferative effects of miR-125a in HCC, promote cell proliferation,35 enhance paclitaxel chemoresistance in HCC cell lines and ia associated with poorer OS.36 These findings suggest that LIN28B may counteract the therapeutic effects of TACE, especially in relation to its chemotherapeutic components, necessitating further research to clarify its role in TACE efficacy. Although these genes are dysregulated in HCC samples, the exact roles and mechanisms by which these genes affect TACE efficacy remain unclear. Further experimental studies are needed to elucidate how their dysregulation influences TACE response and whether targeting these genes could enhance TACE outcomes in HCC.

Pretreatment Images Reveal TFS Gene Expression and TACE Efficacy

Recent studies have explored the predictive potentials of imaging modalities such as contrast-enhanced CT, MRI, and ultrasound-based radiomics to evaluate treatment response and survival in HCC patients undergoing TACE.37–40 While radiomics has shown promise in this area, the biological interpretability of radiomics methodology remains a challenge for its extensive application. In the study by Dai et al,41 radiomics features derived from the arterial and venous phases were employed to construct a model for predicting the efficacy of initial TACE treatment in HCC patients. These radiomic features were found to be related to gene modules that enriched pathways closely linked to tumor development and proliferation, thus establishing a connection between radiomic markers and the biological behavior of tumors. In the present study, the Rad-score was significantly correlated with hub gene expression and TFS score, precisely reflecting TACE refractoriness risk. In addition, external validation results from our cohort demonstrated that the Rad-score could precisely predict responses to initial TACE treatment, with an AUC of 0.676. To further improve the robustness and accuracy of the prediction model for TACE response, baseline and tumor characteristics were also incorporated. The TACE-nonresponders group had greater percentages of patients with a maximal tumor diameter >10 cm and ECOG PS 1 score than did the TACE-responders group. Then, a combined model was established. It demonstrated better performance than the Rad-score or clinical index alone in predicting TACE treatment efficacy. Although Dynamic CT for HCC usually involves the arterial phase, portal vein phase and delayed phase, but only arterial phase images are fully accessible in the TCIA-LIHC database. Despite this limitation, our findings offer another potential rationale for using CT images as a noninvasive tool to predict the efficacy of TACE treatment, especially when combined with biological and clinical markers. However, only CECT images are available from the public database for radiogenomics analysis. Contrast-enhanced computed tomography (CT) and magnetic resonance imaging (MRI) were recommended as the standardized modalities for HCC diagnosis and follow-up. Both modalities have their strengths and limitations. MRI imaging provides superior soft tissue contrast, better functional imaging capabilities without radiation exposure. However, CT imaging is more cost-effective, widely available, suitable for rapid assessment and emergency situations. Further studies about MRI-based radiomics model and MRI versus CT radiomics model for predicting the efficacy of TACE combined with systemic therapies in HCC are needed.

TFS Genes Modulate the Tumor Microenvironment

The tumor microenvironment plays a pivotal role in modulating the development, progression, recurrence and response to immunotherapy in HCC.42 In our study, the high-risk group exhibited an immunosuppressive microenvironment, evidenced by a lower ESTIMATE score, increased infiltration of M0 macrophages, reduced numbers of resting CD4 memory T cells, a lower immune score and an overall immunosuppressive phenotype. The presence of macrophages have been reported to be linked to microvessel density and poor OS and DFS in HCC, and this link is related to their recruitment to tumor tissues and ability to become proangiogenic cells.43 In addition, the density of tumor-infiltrating T cells has been found to be related to prolonged survival in HCC patients,44 which may partially explain the improved OS outcomes in the low-TFS group. This suggests that an active and supportive immune environment, rich in functional T cells, may enhance the efficacy of TACE by promoting immune-mediated tumor clearance. For patients with unresectable HCC who are about to receive systemic treatment, it is difficult to comprehensively determine the characteristics of the tumor microenvironment. If possible, ultrasound- or CT-guided needle biopsy can be considered in routine clinical practice, although these samples do not represent the whole tumor landscape. Further research is needed to investigate the underlying mechanisms and interactions of how hub genes regulate the TME and subsequently induce HCC recurrence or progression after TACE.

TFS Score for Predicting the Response to Targeted Therapy and Immunotherapy

For patients that do not respond to TACE, systemic therapies consisting of tyrosine kinase inhibitors (TKIs), antiangiogenic antibodies, immune checkpoint inhibitors (ICIs), or their combination are recommended.3,45 Several studies have demonstrated that the combination of molecular targeted drugs with ICIs achieves promising objective response rates (ORRs) ranging from 20% to 30% in unresectable HCC patients.46–48 The response of HCC patients to sorafenib can be predicted using our TFS biomarker. Sorafenib responders had higher TFS scores, suggesting that early addition of targeted drugs for patients at high risk of TACE refractoriness might improve their outcomes. Furthermore, lower levels of PD-1 and CTLA-4 expression were observed in the low-risk group, which supported the finding of more favorable immunotherapy responses in the corresponding group. Interestingly, the synergistic effect of TACE with targeted drugs and immunotherapy for HCC has become an area of growing research interest. In advanced HCC, several studies reported that TACE combined with ICIs and anti-VEGF antibody/TKIs as the first-line treatment demonstrated superior tumor control and survival benefits versus TACE or systemic therapies alone.5,6,49 Mechanically, TACE can induce tumor ischemia and hypoxia, which upregulates angiogenic factors like vascular endothelial growth factor (VEGF). Targeted therapies can inhibit VEGF signaling, can prevent tumor revascularization and inhibit further tumor growth after TACE.50 Furthermore, TACE induces necrosis in tumor tissues and causes the release of tumor-associated antigens, then stimulating an immune response, which enhances the activity of ICIs. In addition, TACE modulates the tumor microenvironment by inducing local inflammation, which increases the infiltration of effective T-cells into the tumor, further enhancing the efficacy of immunotherapy.51 In our study, TFS has the potential to serve as an effective indicator for identifying patients likely to respond to sorafenib and immunotherapy. Although the results were based on two another dataset, it helps select individuals benefiting from targeted therapy or immunotherapy. Additionally, several kinds of potential drugs were screened through drug sensitivity analysis in HCC patients with TACE insensitivity, and these drugs may be promising targets for HCC treatment.

Limitations

There are several limitations that should be noted. First, retrospective studies inevitably have inherent bias. In future prospective studies, we will attempt to explore more precise and cost-effective biomarkers based on multiomics integration to aid in the early identification of TACE-insensitive patients. According to the official statement of the GSE104580 dataset, genomic information was obtained from RNA extracted from tumor biopsies of HCC patients treated with TACE and then analyzed using Affymetrix gene arrays. Due to intratumoral heterogeneity, tissues obtained via needle biopsy may not depict the global landscape of HCC lesions. Baseline and clinical information, such as liver function, Barcelona Clinic Liver Cancer (BCLC) stage, presence of portal vein tumor thrombus (PVTT), presence or absence of extrahepatic metastasis, AFP level and hepatitis B virus (HBV) load, and use of molecular targeted drugs or immune checkpoint inhibitors (ICIs), are also important for improving the robustness and prediction performance of the model, but these data are not currently available in public databases. In addition, there are limitations to radiomics studies. Manual delineation of tumors, poor interpretability of radiomics features, and inconsistency of data across different computed tomography (CT) machines may impact the generalizability of the current results. A deep learning-based radiomic model might more intelligently provide clinical guidance and a better bridge between biological features and visible imaging features. Finally, further experiments are required to substantiate the value of the TFS score model.

Conclusion

In conclusion, this study established a gene biomarker capable of predicting the efficacy of TACE for HCC. Additionally, it showed correlations with radiomics features on pretreatment imaging, tumor immune microenvironment characteristics, and the efficacy of systemic therapy in HCC.

Ethical Statement

This study was approved by the Ethics Committee of our center and conducted in accordance with the Declaration of Helsinki. The informed consent was waived due to the nature of retrospective study. All patient records and information were handled with strict confidentiality and privacy. The data used in this study were anonymized to protect patient identities, ensuring that no identifiable information was disclosed at any point in the research process.

Author Contributions

WCD, LB and YR contributed equally to this work. WCD designed the study, performed the experiments, and wrote the manuscript. LB, YR conducted examinations. YZY, LY, DLF, and JH collected the data. CY, YGW and XQY provided hepatological advice and edited the manuscript. CY, XQY revised the manuscript for important intellectual content.

Funding

There is no funding to report.

Disclosure

The authors have no conflicts of interest to declare.

References

1. Brown ZJ, Tsilimigras DI, Ruff SM, et al. Management of hepatocellular carcinoma: a review. JAMA Surg. 2023;158:410–420. doi:10.1001/jamasurg.2022.7989

2. Villanueva A. Hepatocellular carcinoma. N Engl J Med. 2019;380:1450–1462. doi:10.1056/NEJMra1713263

3. Reig M, Forner A, Rimola J, et al. BCLC strategy for prognosis prediction and treatment recommendation: the 2022 update. J Hepatol. 2022;76:681–693. doi:10.1016/j.jhep.2021.11.018

4. Zhai XF, Chen Z, Li B, et al. Traditional herbal medicine in preventing recurrence after resection of small hepatocellular carcinoma: a multicenter randomized controlled trial. J Integr Med. 2013;11:90–100. doi:10.3736/jintegrmed2013021

5. Jin ZC, Chen JJ, Zhu XL, et al. Immune checkpoint inhibitors and anti-vascular endothelial growth factor antibody/tyrosine kinase inhibitors with or without transarterial chemoembolization as first-line treatment for advanced hepatocellular carcinoma (CHANCE2201): a target trial emulation study. EClinicalMedicine. 2024;72:102622. doi:10.1016/j.eclinm.2024.102622

6. Zhu HD, Li HL, Huang MS, et al. Transarterial chemoembolization with PD-(L)1 inhibitors plus molecular targeted therapies for hepatocellular carcinoma (CHANCE001). Signal Transduct Target Ther. 2023;8:58. doi:10.1038/s41392-022-01235-0

7. Kudo M, Matsui O, Izumi N, et al. Transarterial chemoembolization failure/refractoriness: JSH-LCSGJ criteria 2014 update. Oncology. 2014;87(Suppl 1):22–31. doi:10.1159/000368142

8. Tak E, Lee S, Lee J, et al. Human carbonyl reductase 1 upregulated by hypoxia renders resistance to apoptosis in hepatocellular carcinoma cells. J Hepatol. 2011;54:328–339. doi:10.1016/j.jhep.2010.06.045

9. Cheong JH, Wang SC, Park S, et al. Development and validation of a prognostic and predictive 32-gene signature for gastric cancer. Nat Commun. 2022;13:774. doi:10.1038/s41467-022-28437-y

10. He J, Chen Z, Xue Q, et al. Identification of molecular subtypes and a novel prognostic model of diffuse large B-cell lymphoma based on a metabolism-associated gene signature. J Transl Med. 2022;20:186. doi:10.1186/s12967-022-03393-9

11. Wang TC, An TZ, Li JX, Zhang ZS, Xiao YD. Development and validation of a predictive model for early refractoriness of transarterial chemoembolization in patients with hepatocellular carcinoma. Front Mol Biosci. 2021;8:633590. doi:10.3389/fmolb.2021.633590

12. Chen L, Yu CX, Zhong BY, et al. Development of TACE refractoriness scores in hepatocellular carcinoma. Front Mol Biosci. 2021;8:615133. doi:10.3389/fmolb.2021.615133

13. Tang Y, Wu Y, Xue M, Zhu B, Fan W, Li J. A 10-gene signature identified by machine learning for predicting the response to transarterial chemoembolization in patients with hepatocellular carcinoma. J Oncol. 2022;2022:3822773. doi:10.1155/2022/3822773

14. He Q, Yang J, Jin Y. Development and validation of TACE refractoriness-related diagnostic and prognostic scores and characterization of tumor microenvironment infiltration in hepatocellular carcinoma. Front Immunol. 2022;13:869993. doi:10.3389/fimmu.2022.869993

15. Singal AG, Llovet JM, Yarchoan M, et al. AASLD practice guidance on prevention, diagnosis, and treatment of hepatocellular carcinoma. Hepatology. 2023;78:1922–1965.

16. Ligero M, Garcia-Ruiz A, Viaplana C, et al. A CT-based radiomics signature is associated with response to immune checkpoint inhibitors in advanced solid tumors. Radiology. 2021;299:109–119. doi:10.1148/radiol.2021200928

17. Liu Z, Luo C, Chen X, et al. Non-invasive prediction of perineural invasion in intrahepatic cholangiocarcinoma by clinicoradiological features and computed tomography radiomics based on interpretable machine learning: a multicenter cohort study. Int J Surg. 2023. doi:10.1097/JS9.0000000000000881

18. Li Q, Long X, Lin Y, Liang R, Li Y, Ge L. Computed tomography radiomics signature via machine learning predicts RRM2 and overall survival in hepatocellular carcinoma. J Gastrointest Oncol. 2023;14:1462–1477. doi:10.21037/jgo-23-460

19. Wang X, Xie T, Luo J, Zhou Z, Yu X, Guo X. Radiomics predicts the prognosis of patients with locally advanced breast cancer by reflecting the heterogeneity of tumor cells and the tumor microenvironment. Breast Cancer Res. 2022;24:20. doi:10.1186/s13058-022-01516-0

20. Wu J, Liu W, Qiu X, et al. A noninvasive approach to evaluate tumor immune microenvironment and predict outcomes in hepatocellular carcinoma. Phenomics. 2023;3:549–564. doi:10.1007/s43657-023-00136-8

21. Rui W, Zhang S, Shi H, et al. Deep learning-assisted quantitative susceptibility mapping as a tool for grading and molecular subtyping of gliomas. Phenomics. 2023;3:243–254. doi:10.1007/s43657-022-00087-6

22. Kim M, Hui KM, Shi M, Reau N, Aloman C. Differential expression of hepatic cancer stemness and hypoxia markers in residual cancer after locoregional therapies for hepatocellular carcinoma. Hepatol Commun. 2022;6:3247–3259. doi:10.1002/hep4.2079

23. Bossuyt PM, Reitsma JB, Bruns DE, et al. STARD 2015: an updated list of essential items for reporting diagnostic accuracy studies. BMJ. 2015;351:h5527. doi:10.1136/bmj.h5527

24. Thorsson V, Gibbs DL, Brown SD, et al. The immune landscape of cancer. Immunity. 2018;48:812–30.e14. doi:10.1016/j.immuni.2018.03.023

25. Wang H, Shi W, Lu J, et al. HCC: RNA-sequencing in cirrhosis. Biomolecules. 2023;13:141.

26. Liu X, Li T, Kong D, You H, Kong F, Tang R. Prognostic implications of alcohol dehydrogenases in hepatocellular carcinoma. BMC Cancer. 2020;20:1204. doi:10.1186/s12885-020-07689-1

27. Chen Q, Li F, Gao Y, Xu G, Liang L, Xu J. Identification of energy metabolism genes for the prediction of survival in hepatocellular carcinoma. Front Oncol. 2020;10:1210. doi:10.3389/fonc.2020.01210

28. Lau CP, Poon RT, Cheung ST, Yu WC, Fan ST. SPARC and Hevin expression correlate with tumour angiogenesis in hepatocellular carcinoma. J Pathol. 2006;210:459–468. doi:10.1002/path.2068

29. Liu G, Sun J, Yang ZF, et al. Cancer-associated fibroblast-derived CXCL11 modulates hepatocellular carcinoma cell migration and tumor metastasis through the circUBAP2/miR-4756/IFIT1/3 axis. Cell Death Dis. 2021;12:260. doi:10.1038/s41419-021-03545-7

30. Liu G, Yang ZF, Sun J, et al. The LINC00152/miR-205-5p/CXCL11 axis in hepatocellular carcinoma cancer-associated fibroblasts affects cancer cell phenotypes and tumor growth. Cell Oncol. 2022;45:1435–1449. doi:10.1007/s13402-022-00730-4

31. Tokunaga R, Zhang W, Naseem M, et al. CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation—a target for novel cancer therapy. Cancer Treat Rev. 2018;63:40–47. doi:10.1016/j.ctrv.2017.11.007

32. Kinoshita M, Nakamura T, Ihara M, et al. Identification of human endomucin-1 and −2 as membrane-bound O-sialoglycoproteins with anti-adhesive activity. FEBS Lett. 2001;499:121–126. doi:10.1016/S0014-5793(01)02520-0

33. Zhang Q, Wang J, Liu M, et al. Weighted correlation gene network analysis reveals a new stemness index-related survival model for prognostic prediction in hepatocellular carcinoma. Aging. 2020;12:13502–13517. doi:10.18632/aging.103454

34. Zhou X, Liu C, Zeng H, Wu D, Liu L. Identification of a thirteen-gene signature predicting overall survival for hepatocellular carcinoma. Biosci Rep. 2021;41. doi:10.1042/BSR20202870

35. Panella M, Mosca N, Di Palo A, Potenza N, Russo A. Mutual suppression of miR-125a and Lin28b in human hepatocellular carcinoma cells. Biochem Biophys Res Commun. 2018;500:824–827. doi:10.1016/j.bbrc.2018.04.167

36. Tian N, Shangguan W, Zhou Z, Yao Y, Fan C, Cai L. Lin28b is involved in curcumin-reversed paclitaxel chemoresistance and associated with poor prognosis in hepatocellular carcinoma. J Cancer. 2019;10:6074–6087. doi:10.7150/jca.33421

37. Kong C, Zhao Z, Chen W, et al. Prediction of tumor response via a pretreatment MRI radiomics-based nomogram in HCC treated with TACE. Eur Radiol. 2021;31:7500–7511. doi:10.1007/s00330-021-07910-0

38. Sheen H, Kim JS, Lee JK, Choi SY, Baek SY, Kim JY. A radiomics nomogram for predicting transcatheter arterial chemoembolization refractoriness of hepatocellular carcinoma without extrahepatic metastasis or macrovascular invasion. Abdom Radiol. 2021;46:2839–2849. doi:10.1007/s00261-020-02884-x

39. Chen M, Cao J, Hu J, et al. Clinical-radiomic analysis for pretreatment prediction of objective response to first transarterial chemoembolization in hepatocellular carcinoma. Liver Cancer. 2021;10:38–51. doi:10.1159/000512028

40. Liu D, Liu F, Xie X, et al. Accurate prediction of responses to transarterial chemoembolization for patients with hepatocellular carcinoma by using artificial intelligence in contrast-enhanced ultrasound. Eur Radiol. 2020;30:2365–2376. doi:10.1007/s00330-019-06553-6

41. Dai Y, Liu D, Xin Y, et al. Efficacy and interpretability analysis of noninvasive imaging based on computed tomography in patients with hepatocellular carcinoma after initial transarterial chemoembolization. Acad Radiol. 2023;30:S61–S72. doi:10.1016/j.acra.2023.05.027

42. Donne R, Lujambio A. The liver cancer immune microenvironment: therapeutic implications for hepatocellular carcinoma. Hepatology. 2023;77:1773–1796. doi:10.1002/hep.32740

43. Matsubara T, Kanto T, Kuroda S, et al. TIE2-expressing monocytes as a diagnostic marker for hepatocellular carcinoma correlates with angiogenesis. Hepatology. 2013;57:1416–1425. doi:10.1002/hep.25965

44. Garnelo M, Tan A, Her Z, et al. Interaction between tumour-infiltrating B cells and T cells controls the progression of hepatocellular carcinoma. Gut. 2017;66:342–351. doi:10.1136/gutjnl-2015-310814

45. Kudo M, Kawamura Y, Hasegawa K, et al. Management of hepatocellular carcinoma in Japan: JSH consensus statements an