Immune-related cell death index and its application for hepatocellular carcinoma

Identification and validation of immune-related cell death genes in two distinct robust HCC immune subgroups

To investigate the key cell death regulators that participate in immune cell infiltration in TME, we first identified the immune subgroups. Then, the differential expression pattern of these regulators was evaluated in the subgroups (Supplementary Fig. 1). Two to nine immune subgroups were identified by consensus clustering in the TCGA-LIHC cohort based on the relative abundance of 28 defined immune cell types determined by single-sample gene set enrichment analysis (ssGSEA). The consensus matrix, PAC, and cumulative distribution function (CDF) curves confirmed that the optimal number was two (Fig. 1a–d). Meanwhile, the immune subgroups were further confirmed by M3C and NbClust algorithm (Fig. 1b, Supplementary Fig. 2a–c). The two subgroups exhibited substantial differences in immune infiltration, thus regarding as “hot” and “cold” tumors (Fig. 1e and Supplementary Fig. 3a). For validation, the other seven immune estimation algorithms (i.e., ESTIMATE, quanTIseq, xCell, MCP-counter, CIBERSORT, CIBERSORT-ABS and TIMER) were applied, whose results showed the immune infiltration differences between the two subgroups (Supplementary Fig. 3b–h).

To identify the immune-related cell death genes (IRCDGs), we applied differential analysis of all 1078 regulators from 12 diverse cell death modes in the immune-cold- and hot-HCCs (Supplementary Data 1). Among them, 108 differential expression genes (DEGs) were identified and defined as IRCDGs in the TCGA-LIHC cohort (Supplementary Fig. 4a, b, Supplementary Data 3). Principal component analysis (PCA) was performed based on the DEGs expression, and each point was labeled as the identified immune type. Results indicated that the IRCDGs could represent immune subtypes (Fig. 1f). To validate the stability and robustness of the structure of immune subgroups, we performed hierarchical clustering based on the IRCDGs and then obtained two subgroups in each validation cohort. Results showed that the analogous IRCDGs expression patterns and similar proportions of immune subgroups were observed in the validation cohorts (Fig. 1g–i). Submap analysis results also confirmed the similarity of immune subgroups in the validation cohort compared to the referred TCGA-LIHC cohort (Supplementary Fig. 4b).

Compared to previous studies^18,19,20,21, we found that nearly 80% of HCC in both G5 and G6 subtypes from Boyault’s study were classified into the immune-cold subgroups in our research, which could be explained by the over-activation of the WNT signaling pathway^18,22. Meanwhile, G4 subtypes, which had 74% of the immune-hot subgroup and the typical feature of TCF1 mutation, tended to be the hepatocellular adenoma. The other comparison of subtypes with immune subgroups was displayed in Fig. 1j. Overall, HCC immune subgroups in this study were robust and the IRCDGs could represent the immune subgroups tightly associated with immune cell infiltration.

Genomics and functional landscape of immune-related cell death genes and establishment of immune-related cell death index

To further analyze the identified IRCDGs, we performed genomic and function analysis, indicating that IRCDGs consisted of seven cell death modes, including apoptosis, autophagy, necroptosis, etc. (Fig. 2a) The functional annotation of the IRCDGs showed enrichment in the regulation of leukocyte activation, cytokine signaling in the immune system, and regulation of immune effector process pathways (Fig. 2b and Supplementary Fig. 4c). Chromosomal locations and fold-changes for IRCDGs were shown in Supplementary Fig. 4d. Among all the mutations in IRCDGs, missense mutation was the most frequent mutation type. PYD domains-containing protein 3 (NLRP3) and hepatocyte growth factor (HGF) had the highest mutation frequency (6%) (Fig. 2c, d). Based on the copy number variation (CNV) analysis, BIRC3 had the most copy number amplifications, whereas CTSK, ACP5, and CASP4 had the most copy number deletions (Fig. 2e).

a The IRCDGs contain 108 regulators from 7 cell death modes. b Functional pathway network of IRCDGs based on Metascape. c Panel plot displaying frequency summary of genomic mutation types. d Waterfall plot showing the top 20 mutated IRCDGs in HCC patients. e Dumbbell diagram showing gain and loss of CNVs in TCGA-LIHC cohort.

After functional analysis of these IRCDGs, we used a competitive machine learning framework to select the most influential genes together with the respective model to establish the immune cell death index (IRCDI) (Fig. 3a). Of note, nine machine learning models including LASSO (Least absolute shrinkage and selection operator), RSF (Random Survival Forest), Enet (Efficient Neural Network), StepCox, SVM (Support Vector Machine), CoxBoost, SuperPC (Supervised Principal Components), GBM (Gradient Boosting Machine) and plsRcox (Partial Least Squares Regression) were applied and the C-index (Harrell’s concordance index) was calculated to evaluate their performance based on multi-datasets. Results showed although the RSF and GBM had higher C-index in training TCGA cohort, the performance of both decreased in other cohorts (Fig. 3f). Among them, the Cox models, including StepCox and CoxBoost, were more stable and had higher C-index than other models. Finally, the StepCox model (based on nine genes, HSPB8, HMOX1, NQO1, SERPINE1, FGB, FYN, ACSL6, and LGALS9) was confirmed to calculate the immune-related cell death index (IRCDI) (C-index and 95% CI of OS: 0.68 [0.64, 0.73], 0.73 [0.67, 0.79], 0.64 [0.57, 0.71], 0.61 [0.55, 0.67] in TCGA-LIHC, ICGC-LIRI, GSE54236 and GSE14520, respectively). Filtered genes by univariate Cox analysis were shown in Fig. 3b. And Fig. 3c–e showed the LASSO regression progression and coefficients of these genes. In the TCGA-LIHC cohort, significant stratification of groups was achieved on all model genes except the LGALS9 by the Log-Rank test, and the KM curves were shown in Supplementary Fig. 5. For benchmark survival analysis, the HRs of these model genes were calculated based on nine independent datasets and the results were shown in Supplementary Fig 6a–i. The heatmap in Fig. 3h visually depicted the scaled expression patterns of the model genes and clinical features in each HCC patient. Furthermore, we developed a free online website ( for the easy application of the IRCDI model, which allowed clinicians to input the value of model genes and quickly obtain the predicted results (Fig. 3g).

a Competitive machine learning framework for developing the IRCDI. b HRs with 95% CIs of 29 survival-related genes by univariate Cox analysis. c, d Determination of optimal variables based on minimum lambda. e The coefficients of nine model genes were identified using the StepCox model. f A combined heatmap shows the C-index of each machine learning method in each dataset. g Screenshot of the developed IRCDI online website. h Heatmap showing the relationship between IRCDI and clinical features of model genes.

Clinical associations and survival prediction performance of IRCDI

After constructing the IRCDI, we initially investigated the correlation between clinical features and IRCDI using the TCGA-LIHC cohort. Our results demonstrated that IRCDI was not influenced by baseline clinical features such as gender and age (Fig. 4a, b). Interestingly, high-grade HCC exhibited significantly higher IRCDI values compared to low-grade HCC (G1 vs G3, p = 0.0022; G2 vs G3, p = 0.0038; Fig. 4c). Additionally, HCC with high-IRCDI displayed a greater extent of residual tumor and elevated levels of alpha-fetoprotein (Fig. 4d, e). Notably, IRCDI exhibited strong associations with tumor malignancy and prognosis. Patients who succumbed during the follow-up period exhibited higher IRCDI values compared to those who survived (Fig. 4f). Moreover, patients in the T2 and T3 stages exhibited higher IRCDI values than those in the T1 stage. These findings were further supported by the American Joint Committee on Cancer (AJCC) stage classification (Fig. 4g, h). Furthermore, the relationships between IRCDI and pN, and pM were also investigated (Fig. 4i, j).

a–j Boxplots for gender, age, grade, residue tumor, AFP, OS, stage, pT, pN, and pM, respectively. Centre line: the median; Bounds of the box: the 25th and 75th percentile; k Survival curves of IRCDI of TCGA-LIHC OS (Left) and forest plot of multivariate Cox regression results (Right). l Survival curves of IRCDI of TCGA-LIHC RFS (Left) and forest plot of multivariate Cox regression results (Right). m–p Survival curves of IRCDI of LIGC-LIRI OS, GSE54236 OS, GSE14520 OS, and GSE14520 RFS, respectively. q–t Survival curves of IRCDI of TCGA-PAAD OS, TCGA-STAD OS, TCGA-ESCA OS, and TCGA-COAD OS, respectively. Statistical significance: *P 

To assess the independent prognostic value of IRCDI in HCC, we conducted multivariate Cox regression analysis, incorporating important clinical features. Our findings revealed that IRCDI remained a significant predictor for overall survival (OS) and recurrence-free survival (RFS) (Fig. 4k, l). To evaluate the predictive accuracy of 1-, 3-, and 5-year survival, we employed Kaplan–Meier survival analysis. Notably, all p-values from the log-rank test were less than 0.001. In the TCGA-LIHC cohort, the AUCs with 95% CIs for 1-, 3-, and 5-year OS were 0.73 [0.65, 0.81], 0.71 [0.66, 0.80], and 0.68 [0.60, 0.76], respectively (Supplementary Fig. 7a). Furthermore, the ICGC-LIRI and GSE54236 cohorts demonstrated improved 1- and 3-year OS predictions (ICGC-LIRI: 1-year OS AUC, 0.78 [0.66, 0.80]; 3-year OS AUC, 0.7 [0.60, 0.80]; GSE54236: 1-year OS AUC, 0.79) (Fig. 4m, n and Supplementary Fig. 7c, d). These results validated the robust performance of the IRCDI model in predicting OS. Although the model’s performance in RFS prediction was slightly lower than that of OS, it still exhibited predictive capability (TCGA-LIHC: 1-year RFS AUC, 0.68 [0.61, 0.74]; 3-year RFS AUC, 0.64 [0.56, 0.72]; 5-year RFS AUC, 0.67 [0.57, 0.77]). Comparatively, our model showed better performance than the clinical stage and other features. The calibration curve and C-index of the other models can be found in Supplementary Fig. 7g, h, respectively. To fully validate the applicability of the IRCDI, we also assessed its performance in other digestive tract tumors. Our results demonstrated a significant survival difference between high and low-IRCDI groups in all cohorts, including TCGA-PAAD, TCGA-STAD, TCGA-ESCA, and TCGA-COAD (Fig. 4q–t). Since the clinical values of IRCDI had been illustrated, we then explored the distribution of each model gene in the TNM stage, T, N, and M stage in the TCGA-LIHC cohort. The results are displayed in Supplementary Fig. 8a–i. In conclusion, our developed IRCDI model is user-friendly and has the potential for widespread application among clinicians.

Different driver genes and activated biological pathways of IRCDI subgroups

To investigate the related biological processing of IRCDI, we first divided patients into two subgroups according to the median IRCDI in the TCGA-LIHC cohort. Of these, 42% of high-IRCDI patients had alterations on the driver gene TP53, whereas the most frequently altered gene was CTNNB1 (36%) in low-IRCDI patients (Fig. 5a, b). Statistical significance of these driver genes was confirmed, which indicated the carcinogenic mechanisms (Fig. 5c, e and Supplementary Fig. 9a). Meanwhile, both TP53 and CTNNB1 mutation combined with the IRCDI showed different risk layers (Fig. 5d, f). We then explored the expression of each gene in TP53 cancer pathways. The results showed that the TP53 cancer pathway was more active in patients with high-IRCDI with the positive Spearman correlation (Supplementary Fig. 9c, d).

a, b Top 10 mutated genes in IRCDI subgroups based on TCGA-LIHC cohort. c Lollipop Plot illustrating the TP53 gene mutated sites and types in IRCDI subgroups. d Survival curve demonstrating different risk layers based on IRCDI and TP53 mutation status. e Lollipop Plot showing the CTNNB1 gene mutated sites and types in IRCDI subgroups. f Survival curve demonstrating different risk layers based on IRCDI and CTNNB1 mutation status. g Gene Set Enrichment Analysis (GSEA) showing the up- and down-regulated pathways in IRCDI subgroups. h Scatter plot displaying the most significant active pathways related to IRCDI. i, j Box plots showing the differences in representative pathways. Centre line: the median; Bounds of the box: the 25th and 75th percentile; k The relationship between module genes and IRCDI and clinical features based on WCGNA analysis. l Pathway enrichment analysis of turquoise and green module genes.

Then, we analyzed the cancer patterns in IRCDI groups. Gene set enrichment analysis (GSEA) analysis indicated that compared to the low-IRCDI HCC, the high-IRCDI HCC had the activated cell cycle pathways and suppressed metabolism pathways (Fig. 5g). The activated cancer pathways were significantly different in distinct IRCDI groups (Supplementary Fig. 9b). The spearman correlation analysis revealed that, in consensus with the GSEA results, cell cycle-related pathways including the well-known MYC/E2F target, G2M checkpoint and mitotic spindle pathways were significantly activated in high-IRCDI group (Fig. 5h, i). Meanwhile, dysregulated metabolism pathways, including bile acid, xenobiotic, fatty acid, and heme metabolism, were observed in the low-IRCDI group (All p-values 5j). In addition, the cell proliferation pathways with the key gene CDK1 were significant activated in IRCDI subgroups with good discriminability (Supplementary Fig. 9e, f). For support, a widely used marker used in clinical, MIK67, is analyzed, showing higher expression in the high-IRCDI group with an AUC of 0.75 (Supplementary Fig. 9g, h). In order to identify the most significant gene cluster with IRCDI, we performed a weighted correlation network analysis (WGCNA) analysis in the TCGA-LIHC cohort (Supplementary Fig. 10a–c). Ten gene modules were clustered to investigate the correlation between IRCDI and clinical features (Fig. 5k). Results indicated that the turquoise module had a significantly positive correlation to IRCDI and AJCC stage (all p-values 10d, f, respectively. Functional enrichment results showed turquoise module genes were enriched in cell cycle and DNA replication pathways and green module genes were enriched in the metabolism of amino acids and derivatives and small molecule catabolic process pathways. Other enrichment results were displayed in Fig. 5l and Supplementary Fig. 10e, g. The above results explained that patients with higher IRCDI tended to have more activated cell proliferation pathways, which could lead to a worse prognosis.

IRCDI is related to immune cell infiltration and PD-1 expression level

To investigate the IRCDI distribution in different cell types of TME, we performed single-cell analysis from the primary HCC without any cell type enrichment. By the singleR package, cells were divided into five clusters (Fig. 6a). The IRCDIs of each cluster were calculated and displayed in Fig. 6b. Results showed that myeloid cells had the highest IRCDI, followed by fibroblast cells (Fig. 6c). The results provided biological evidence for the application of the model to predict the prognosis of HCC patients.

a UMAP analysis showing identified cell types. b, c Distribution of IRCDI across different cell types. d Boxplot demonstrating the differential IRCDI levels in experimental immune subtypes. Centre line: the median; Bounds of the box: the 25th and 75th percentile; e Immune cell infiltration in IRCDI subgroups. Statistical significance: *P f–h Scatter plot illustrating the negative correlations between activated CD8 T cells, effector memory CD8 T cells, and type 1T helper cells. i Bar plot depicting the relationship between CD8A and IRCDI based on proteomics. j–m Boxplot demonstrating the negative correlation between PD-1 expression and IRCDI in TCGA-LIHC, ICGC-LIRI, GSE14520, and GSE54236 datasets. Centre line: the median; Bounds of the box: the 25th and 75th percentile; n Immunohistochemistry experiment validating the association of PD-1 expression with IRCDI.

In this study, HCC immune subtypes were defined as immune-hot and immune-cold, which are usually called immune-desert and immune-inflamed. In fact, another immune subtype named immune-excluded was well-known, which has immune cell infiltration but low immune response because of the immune checkpoint effects. We used the IMvigor 210 cohort to investigate whether the IRCDI could distinguish the above immune subgroups. Results showed that patients in different immune subgroups gained different IRCDI scores. The immune-desert group had a higher IRCDI than the immune-excluded group and immune-inflamed group (p = 0.045 and 0.00055, respectively; Fig. 6d). Out of expected, even though we didn’t use any specific information about the immune-excluded subtype for training the model, the IRCDI model still performed well in distinguishing the immune-excluded patients from the immune-inflamed patients (AUC: 0.6024, Supplementary Fig. 11c). We then investigated the correlations between IRCDI and 28 kinds of immune cells. Many immune cells infiltered varied in high- and low- IRCDI groups (Fig. 6e). The correlation analysis results showed activated CD8+T cell, effector memory CD8+T cell, and type 1T helper cell were significantly negatively with IRCDI (r = −0.2008, −0.2021 and −0.3052; p = 2e-04, 2e-04, 0; Fig. 6d–h). Other immune cell results were displayed in Supplementary Fig. 11a. For validation, although the statistical result was not significant because of the small samples, the protein level of CD8A was negative correlated with the IRCDI in our cohort (r = −0.4286, p = 0.3536; Fig. 6i).

We further explored the differential activation of each immune step in low- and high-IRCDI subgroups. Results indicated that the high-IRCDI group tended to recruit more MDSC cells, which had been confirmed as functioning in the immune response suppression (p 12)²³. Results also indicated that there was no difference in effector T cells recruiting, while the processing of infiltration of immune cells into tumors varied in IRCDI subgroups (p 12). Apart from immune cell infiltration, another factor that determined the immune inhibitors efficacy was the immune checkpoint expression level. Notably, PD-1 showed significantly low expression in the high-IRCDI group in independent datasets (Fig. 6j–m; p = 0.015, 0.034, 0.057, and 0.016 in TCGA-LIHC, ICGC-LIRI, GSE14520, and GSE54236, respectively). The results were also validated by IHC in our cohort (Fig. 6n). All the above results explained that the patients with higher IRCDI tended to have a less immune cell infiltrated TME.

Application of IRCDI to predict immunotherapy and trans-arterial chemoembolization response

Several treatments have been proposed for advanced hepatocellular carcinoma (HCC), including immunotherapy, targeted drugs, and trans-arterial chemoembolization (TACE). These treatments have been shown to induce tumor cell death, which was subsequently cleared by immune cells²⁴. Because the IRCDI had strong associations with the tumor cell death progression and immune microenvironment, we then applied TIDE analysis to validate IRCDI, four real-world immunotherapy cohorts and one TCAE treatment cohort for a full evaluation of the IRCDI for drug response prediction (Fig. 7a). Our results showed that TIDE scores (tumor immune escape score) increased with IRCDI, and were significantly different in high- and low- IRCDI groups (Fig. 7b). Survival analysis also indicated that patients with high-IRCDI under immune checkpoint inhibitor treatment had the worse prognosis (Log-rank p-values: 0.0003 and 0.049; Fig. 7c, d). Then, we compared IRCDI to the well-known marker PD-1 in the GSE35640 and GSE91061 cohorts receiving anti-MAGE-A3 and anti-PD-1/CTLA-4 immunotherapy to estimate the prediction accuracy of the IRCDI model. Results showed that IRCDI had a better performance than PD-1 in GSE35640 cohort (AUC: 0.76 vs 0.51, Fig. 7e). Another result indicated that IRCDI had nearly the same performance compared to PD-1 in GSE9106 (AUC: 0.65 vs 0.69, Fig. 7f). In the TACE treatment cohort, IRCDI still performed well in predicting the response, while PD-1 performed very poorly (AUC: 0.75 vs 0.51). These results indicated that most typical markers only performed well under specific conditions. Compared to the typical markers, IRCDI demonstrated strong generalization capability and was more suitable for clinical application.

a Flow chart of treatment prediction. b Boxplot demonstrating the differential TIDE scores in IRCDI subgroups (Left). Scatter plot showing the linear relationship between the IRCDI and TIDE score (Right). c, d Overall survival curves of IRCDI under the ICB treatment in GSE100797 and IMvigor 210 cohort. e Box plots and ROC curves show the performance of anti-MAGE-A3 response prediction in GSE35640 of IRCDI (Left panel) and PD-1 (Right panel), respectively. f Box plots and ROC curves show the performance of anti-PD-1/CTLA-4 response prediction in GSE91061 of IRCDI (Left panel) and PD-1 (Right panel), respectively. g Box plots and ROC curves show the performance of TACE response prediction in GSE10458 of IRCDI (Left panel) and PD-1 (Right panel), respectively. Centre line: the median; Bounds of the box: the 25th and 75th percentile.

Exploration of available drugs according to IRCDI

Next, we screened the sensitivities of the common anti-cancer drugs according to the calculated IRCDI to rescue the situation of ICI-resistant patients. A flow chart is displayed in Fig. 8a. The predicted IC50 values of nearly two hundred FDA-approved anti-cancer drugs were calculated using the “oncoPredict” package. Resistance and sensitivity of typical drugs in patients with high-IRCDI are shown in Fig. 8b. Then, we calculated the IRCDI in each HCC cell line and used the GDSC (Genomics of Drug Sensitivity in Cancer) database to validate the predicted drugs. The intersected drugs were regarded as the available drugs, together with whose targets were displayed in Fig. 9c. The involved pathways were shown in Fig. 8d (Gene count >1). Among the five available drugs, sepantronium bromide and AZD7762 were selected for further analysis. The two drugs not only had different IC50 in the IRCDI HCC cohort but also had significant negative correlations in HCC lines confirmed by lab experiments (Fig. 8e–g, and Supplementary Fig. 13a–c). Further analysis showed the drug effects depended on the expression level of their target BIRC5 and CHK1(Fig. 8i and Supplementary Fig. 13e, f). Both the targets were associated with HCC prognosis (Log-rank p-values: 6.7e-05 and 3.1e-105; Fig. 8h and Supplementary Fig. 13d). According to the value of IC50, sepantronium bromide was regarded as the most potential drug. We then used Autodock to calculate the binding site of the drug and the identified target. Results of the exact binding site and binding energy (−4.409 Kcal/mol) were displayed in Fig. 8i.

a Flow chart of drug selection. b Spearman analysis showing the top 10 resistance and sensitivity drugs of IRCDI. c Venn diagram displaying the overlap drugs of the predicting results and GDSC recorded. d Venn diagram displaying the overlap drugs of the predicting results and GDSC recorded. e The IC50 values of sepantronium bromide in each HCC cell line were ordered by IRCDI. f, g Boxplot showing BIRC5 expression level and predicted IC50 of sepantronium bromide in IRCDI subgroups. Centre line: the median; Bounds of the box: the 25th and 75th percentile; h Survival curve of BIRC5 in TCGA-LIHC cohort. i Scatter plot showing the relationship between the IRCDI and the predicted IC50 of the drug and the target gene, respectively. j The 3D structure shows the binding pattern of the selected drug and target according to the Autodock results.

a Flowchart of the validation experiments. b QRT-PCR results showing the ACSL6 differential expressed in normal and cancer tissue. c IHC results from the HPA database. d Bar plot showing the relative ACSL6 mRNA level in control and Si-ACSL6 HepG2 cell line. e Cell proliferation curves of control and Si-ACSL6 HepG2 cell line. f Bar plot showing the relative ACSL6 mRNA level in control and Si-ACSL6 Huh7 cell line. g Cell proliferation curves of control and Si-ACSL6 Huh7 cell line. h, i Wound healing assay results of control and Si-ACSL6 HepG2 cell line. j, k Wound healing assay results of control and Si-ACSL6 Huh7 cell line. l Pathway enrichment results of different expressed genes by proteomics of control and Si-ACSL6 HepG2 cell line. m–p RFTN1, URGCP, LETMD1, L1CAM protein level in control and Si-ACSL6 HepG2 cell line. Statistical significance: *P 

Experimental validation of IRCDI panel genes by in-house samples and HCC cell lines

The prognosis and drug response prediction performance of IRCDI had been illustrated, and for the further investigation of the model genes, we collected seven HCC tumors and normal samples and applied qRT-PCR and LC–MS/MS experiments. Especially we validated the key model genes by siRNA and proteomics assay (Fig. 9a). As shown in Fig. 9b and Supplementary Fig 14a–h, ACSL6, FGB, and FYN were significantly downregulated in cancer samples, while HSPB8 and NQO1were up-regulated. Of which, highly expressed ACSL6 tended to have a lower stage and pT (Supplementary Fig. 8a). For support, the IHC result confirmed that ACSL6 expressed lower in cancer compared to the normal tissue (Fig. 9c). Our research findings strongly supported the negative correlation between low ACSL6 expression and prognosis in patients diagnosed with HCC. We then calculated cancer pathways in our proteomics data. After that, spearman correlation analysis was used to evaluate the relationship between the expression of ACSL6 and cancer pathways. Results indicated that most cancer pathways were negatively correlated with the expression of ACSL6, especially the cancer metastasis and stemness pathways (r: −0.8768 and −0.9536; p: 0.0095 and 0.0008, respectively. Supplementary Fig. 14i–k).

To further validate the impact of ACSL6 on HCC, we conducted additional in vitro experiments. Using siRNA-mediated knockdown, we successfully reduced the mRNA levels of ACSL6 in HepG2 and Huh7 cells (Fig. 9d, f). The CCK-8 cell proliferation assays demonstrated a significant enhancement in cell growth for both HepG2 and Huh7 cell lines upon ACSL6 reduction (Fig. 9e, g). Similarly, the downregulation of ACSL6 resulted in augmented migration of HCC cells, confirmed by scratch wound healing assay (Fig. 9h–k). The proteomics results revealed that compared to the control groups, si-ACSL6 groups had the deregulated protein localization and autophagy biological processing (Fig. 9l). Of note, genes of NF-kB and EMT pathways, including RFTN1, LETMD1, URGCP, and L1CAM in Si-ACSL6 groups were up-regulated (Fig. 9m–p). These results provided compelling evidence supporting the crucial role of ACSL6 as a promising candidate for a prognostic biomarker in HCC.