Construction of an Immunogenic Cell Death-Related Gene Signature and Genetic Subtypes for Predicting Prognosis, Immune Microenvironments, and Drug Sensitivity in Hepatocellular Carcinoma

Introduction

The prevalence of liver cancer is rising globally, and it continues to be a problem for global health.^1,2 By 2025, liver cancer is expected to impact over one million people.³ Hepatocellular carcinoma (HCC) is the most common type of liver cancer, accounting for 90% of cases.⁴ The development of HCC is a complex process involving ongoing inflammatory damage, regeneration, and necrosis of hepatocytes. The survival of cancer cells was facilitated by anomalies in molecular signaling and gene expression, which was the core reason for the development of cancer.⁵

In recent years, immunotherapy based on the modulation of the tumor immune microenvironments has become a new therapy for HCC patients, but remission rates were unsatisfactory. HCC was a typical immunogenic cancer that usually develops in a setting of ongoing inflammation.⁶ HCC developed while the immune system was in a “cool” state, shielding it from cytolytic attack by lymphocytes entering the tumor. The purpose of immunotherapy was to transform the “cold” state of tumor immunity into a “hot” state and thus restructure the tumor immune environment to enhance the anti-tumor immune responses.⁷ The main difficulty with immunotherapy for HCC was that the immune system was prevented from exerting its typical anti-tumor effects by the tumor immune escape mechanisms, the tumor suppressive immunological microenvironments and the low immunogenicity of tumor cells. Exploring the precise biomarkers for immunotherapy will be especially crucial given the advancements in immunotherapy.

Immunogenic cell death (ICD) was thought to be one of the most effective methods for completely eliminating tumor cells due to its ability to stimulate immune responses and create long-lasting immunological memory.^8–10 Increasing evidence demonstrated that some chemotherapy drugs not only directly killed tumor cells, but also had the capacity to induce ICD.^11–14 ICD was a type of cell death caused by releasing tumor associated antigen (TAA) and tumor specific antigen (TSA) in order to activate the immune response of the immune system.¹⁵ The immunogenic characteristics of ICD were mainly mediated by damage-associated molecular patterns (DAMPs), which include surface-exposed calreticulin, secreted ATP and released high mobility group protein B1.⁹ The production of ICD-related DAMPs and other immunostimulatory molecules facilitated the formation of efficient connections between innate immune cells (such as DCs or macrophages) and cancer cells, thereby initiating a treatment-related adaptive immune response.¹⁶ Due to the characteristics of immunogenicity production and tumor antigen presentation, ICD-related genes (ICDRGs) were expected to provide new perspectives and methods for anti-tumor immunotherapy. Therefore, identifying the signature of ICDRGs in HCC is essential for the prognosis and treatment, which allows us to understand the immune microenvironments better and utilize immunotherapy of HCC.

Bioinformatics tools and techniques enhance the precision and efficiency of transcriptomic feature analysis, enabling comprehensive comprehension of gene expression patterns and regulatory mechanisms.¹⁷ Transcriptomic feature analysis was utilized to categorize cancer into clinically significant molecular subtypes for diagnosis, prognosis, or therapeutic selection purposes.¹⁸ By analyzing gene expression patterns in tumor tissues or cells, distinct subtypes were identified, which may exhibit significant differences in cancer progression, survival rates, and treatment responses. This approach not only provides clinicians with more accurate survival prediction information but also facilitates the development of personalized treatment strategies, thereby improving treatment efficacy and patient survival rates. One of the primary methods for predicting the survival prognosis in cancer patients is to construct risk prediction models based on key genes.^19–21 The combination of both not only provides more accurate survival prediction information for clinical doctors but also assists in formulating individualized treatment strategies, thereby improving patient treatment outcomes and survival rates.

In this study, we constructed genetic subtypes and risk models based on ICDRG expression characteristics to further investigate the tumor immune landscape, prognosis, and clinical benefit of immunotherapy and chemotherapy in HCC patients. These results will provide a comprehensive understanding of the impact of ICDRGs in the development of HCC and help improve the efficacy of individualized treatment for HCC patients.

Materials and Methods

Data Collection and Gene Identification

The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/repository) was used to acquire RNA sequencing (RNA-seq) expression data, genomic mutation data, and clinical data for HCC patients. The International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/projects/LIRI-JP) was also used to retrieve the RNAseq expression data and corresponding clinicopathologic data. Fifty-five ICDRGs were collected from extensive previous studies^22,23 (Table S1). The analysis of the protein-protein interactions (PPI) was carried out by STRING (https://string-db.org). Figure 1 illustrated the flowchart for the study design.

Figure 1 The flowchart of this study.

Consensus Clustering Analysis

The ConsensusClusterPlus algorithm²⁴ was used to reclassify patients in order to explore the heterogeneity of ICDRG expression in HCC. For the purpose of confirming the ideal number of subtypes, the cumulative distribution function (CDF) and consensus matrices were used.

Tumor Mutation Burden Analysis

The mutation landscape was plotted as a waterfall using the R package “maftools”²⁵ showing the genes with the highest mutation frequency (Top 20). The copy number variation (CNV) frequencies were calculated, and the above results were shown in lollipop plots. The locations of these genes on the chromosomes were visualized using the “RCircos” package of R software.

Immunogenomic Landscape Analysis

The “estimate” package utilized the ESTIMATE algorithm²⁶ to calculate the immunological or stromal fraction of the tumor microenvironments (TME). In order to evaluate the immune-cell invasion, the CIBERSORT algorithm (http://cibersort.stanford.edu/) was employed. The ssGSEA analysis was used to compute ssGSEA scores that reflected the frequency of each immune cell type.²⁷ Immune checkpoint and human leukocyte antigen (HLA) gene expression levels were compared in the different groups.

Differential Gene Expression (DEG) Analysis and Functional Pathway Enrichment Analysis

DEG analysis was carried out using Limma (cutoff of FDR ≤ 0.05; fold change ≥ 1.5). Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analysis (GSEA) was used to assess the biological roles of the prognostic genes by the “clusterProfiler” R package.²⁸

Construction and Evaluation of the Risk Model

Differentially expressed genes between tumor and normal tissue were identified using the RNA-seq results. We developed a Cox regression model with Least Absolute Shrinkage and Selection Operator (LASSO) algorithms feature selection to analyze the roles of ICDRGs in survival prediction. The TCGA and ICGC cohorts were used, respectively, as training sets and validation sets. The discrimination and calibration of the risk model were assessed using receiver operating characteristics (ROC) curves and calibration curves. To assess the model’s diagnostic value and applicability, clinical impact curves (CIC) were performed by using the resample bootstrap method (No. of bootstrap replications = 1000). The R package “regplot” was used to create a predictive nomogram combining clinical variables (age, gender, stage, and grade) and ICDRGs risk. Uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding (t-SNE) were used to perform a dimensional reduction analysis and visualization. The risk score and vital status of each person were displayed on the risk curves and scatter plots. The differentiation in OS between patients at low- and high-risk were analyzed using Kaplan-Meier (K-M) survival curves. Sankey plots illustrated the connections between patients’ risk scores, ICDRG clusters, and survival status.

Cell Culture

HCC cells (Huh7 and HepG2) and normal adult liver epithelial THLE-2 (as control cells) were obtained from Fubo Bio (Beijing, China) and maintained in Dulbecco’s modified eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100 μg/mL streptomycin. The 5×10⁵ LO2, HepG2, and Huh7 cells were seeded in 6-well plates at 37°C in 5% CO₂ with saturated humidity.

Immunohistochemistry Analysis

The expression pattern of the ICDRGs at the protein level was examined using the Human Protein Atlas database (HPA, http://www.proteinatlas.org), which contained protein data from tumor and normal clinical samples. Pathology and tissue sections, respectively, provided photomicrographs of ICDRGs immunohistochemistry (IHC) staining in HCC and matching normal tissues.

Chemotherapeutic Drug Sensitivity Analysis

We assessed the 50% inhibitory concentration (IC50) of medications using the pRophetic algorithm to ascertain the relevance of the risk score in forecasting chemotherapies and targeted drug sensitivity. The cgp2016 dataset used for prediction included 251 drugs taken from the Genomics Database for Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/). The R package “pRRophetic” was used to calculate IC50 by using ridge regression. The default settings were used for each parameter. Using the R package “ggpubr”, we compared the difference in IC50, and the results were shown as boxplots.

Statistical Analysis

The statistical analysis for this study was performed using R software (version 4.2.0), and the bioinformatics analysis was facilitated by the Sangerbox platform.²⁹ The Wilcoxon rank-sum test was used to compare variables that were not normally distributed, and the independent Student’s t-test was used to compare continuous variables between two groups. Categorical variable data were compared using the chi-squared test. Correlations were analyzed using the Pearson chi-square test. Differences were considered statistically significant when the P-value was < 0.05.

Results

Identification of ICDRGs and Four Subtypes in HCC Patients

Protein-protein interaction (PPI) analysis further revealed the association between 55 ICDRGs (Figure 2A). 35 ICDRGs showed different expression patterns between normal and HCC patients, of which 19 were upregulated and 16 were downregulated in HCC, respectively (P < 0.05, Figure 2B). The k-means consensus clustering analysis was used to classify the TCGA cohort into four different subtypes (Figure 2C–E). As shown in Figure 2F and G, we found that the vast majority of ICDRGs were upregulated in cluster 3(C3), corresponding to its median survival time being the shortest among the four groups. Cluster 2 (C2), on the other hand, had the longest median survival time, followed by Cluster 1 (C1), Cluster 4 (C4), and C3 (P = 0.022, K-M survival analysis), while the majority of the ICDRGs in C2 were down-regulated. The results suggested that patients with low expression levels of ICDRGs had prolonged median survival time, whereas those with high expression levels of ICDRGs had shorter median survival time, suggesting a poor prognosis.

Figure 2 Identification and prognostic evaluation of four subtypes based on ICDRGs in HCC. (A) The PPI network of ICDRGs. (B) The heatmap of differential gene expression between tumor and normal tissues. (C) Consensus clustering cumulative distribution function (CDF) for k = 2 to 9. (D) Relative change in the area under the CDF curve for k = 2 to 9. (E) HCC patients in the TCGA cohort were divided into four distinct clusters when k = 4. (F) The expression of ICDRGs expression in four subtypes (G) K-M survival analysis of the OS status of HCC patients in four ICDRGs subtypes.

Tumor Mutation Burden and Tumor Immune Microenvironments in Four ICDRG Subtypes

The mutation frequencies of the ICDRGs (top 20) in HCC were initially identified in the four clusters (Figure 3A–D). The mutation analysis revealed four subtypes with varying characteristics and degrees of mutations in HCC patients, which were mainly concentrated in four genes: TP53, CTNNBA, TTN, and MUC16. The mutation of CTNNB1 was the most important mutation event in C2, accounting for 32% of the total mutation events. However, the CTNNB1 mutation only accounted for 12% of the mutational events in C3, which had the lowest percentage of such mutation events among the four groups. Comparatively, the frequencies of TP53 mutations were lowest in C2 (18%) and highest in C3 (34%). The majority of gain-of-function mutations were located in TLR3, TP53, HMGB1, and ATG5, while the majority of loss-of-function mutations were located in NLRP3, IL10, PIK3CA, and TLR5 (Figure 3E–H). The locations of these genes on the chromosomes were shown in Figure 3I–L. Then, we explored the heterogeneity of the immune microenvironments in the different subtypes. Overall, C3, the cluster with higher expression of ICDRGs, showed lower tumor purity, higher stromal scores and higher immune scores compared to other subgroups (Figure 4A and B). In addition, most immune checkpoints were significantly upregulated in C3 (Figure 4C). The CIBERSORT algorithm was used to compare the immunological infiltration of 22 immune cells among the four subtypes (Figure 4D and E). There were significant differences in the percentage of plasma cells, CD4⁺ T cells memory resting, CD4⁺ T cells memory activated, follicular helper T cells, regulatory T cells (Tregs), NK cells activated, Macrophages M0, Mast cells and Neutrophils (P < 0.05). The subtype C3, which exhibited higher expression of ICDRGs, produces more immune cells that produce antibodies or mediate related responses (plasma cells and follicular helper T cells). However, the immune cells with phagocytic and cytotoxic functions were reduced in these patients. There was an increase in the expression of macrophage M0 and a decrease in the expression of macrophage M1, M2 and activated NK cells. The subtype C2 as a cluster with lower ICDRGs expression showed the opposite results to C3. These results suggested that the expression of ICDRGs might have affected the prognosis by influencing the number of immune cells.

Figure 3 Tumor mutation burden of four ICDRG subtypes in HCC patients. (A-D) waterfall maps of the somatic mutation landscape in four subtypes. (E-H) lollipop diagrams of the copy number abnormalities indicated the degree of copy number loss (green) or gain (red). (I-L) Circus plots of the chromosome distributions of 35 differentially expressed genes. The outermost ring represented the 23 human chromosomes, while the inner ring displayed the 35 genes linked to different chromosomes. Red, blue, and black dots indicated copy number gain, loss, and no change, respectively.

Figure 4 Immune microenvironment of four ICDRG subtypes in HCC patients. The stromal scores, immune scores (A), and tumor purity (B) for four subtypes by Mann–Whitney U-test. (C) The boxplot showed the expression of immune checkpoints in four subtypes. The boxplot (D) and heatmap (E) of immune infiltration cells between four subtypes of HCC. (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

Establishment and Validation of the ICDRG Risk Models in HCC

Univariate Cox regression analysis revealed that eleven ICDRGs were associated with overall survival (OS) (P < 0.05, Figure 5A). We created a prognostic model based on ICDRGs. LASSO Cox regression was used to further filter the genes and reduce the risk of overfitting (Figure 5B and C). Eight genes (ATG5, BAX, CASP8, HSP90AA1, LY96, NT5E, P2RY2 and TLR2) were identified for the construction of the ICDRG risk model, and a nomogram with the risk scores, clinical features, and pathological indicators was created (Figure 5D). Next, We performed a time-dependent ROC analysis. In the TCGA cohort, the AUC was 0.731, 0.657, and 0.666 for the 1-year, 3-year and 5-year OS, respectively (Figure 5E). In the ICGC cohort, the AUC values were 0.626, 0.630 and 0.686 for the 1-year, 3-year and 5-year OS, respectively (Figure 5F). The AUC of the ICDRG risk scores in the TCGA cohort was significantly higher than age, gender, tumor stage, and pathological grade, demonstrating better discrimination of the ICDRG risk model (Figure 5G and H). UMAP and t-SNE analysis indicated that two independent subgroups of HCC were well-distinguished using the ICDRG risk model (Figure S1). The calibration curves showed good agreement between the predictions of the model and the actual observed probabilities (Figure 5I). The results of the clinical impact curves showed good performance for the model according to clinical application (Figure 5J and K). As the risk threshold increased, there was a decrease in unnecessary treatment and an increase in net clinical benefit. Then, we calculated scores for each patient using the ICDRG risk model. The risk score of each HCC patient was calculated as follows: risk score = 0.1946*ATG5 + 0.1023*BAX + 0.0394*CASP8 + 0.2425*HSP90AA1 + 0.0320*LY96 + 0.0135*NT5E + 0.2086*P2RY2 + 0.0138*TLR2. As demonstrated in Figure 5L and M, univariate and multivariate Cox regression analysis revealed that the risk scores were independent predictors of OS compared to other clinical indications. The results suggested that the ICD risk scores hold significant prognostic value in predicting the survival outcome of the patients, independent of other clinical factors such as age, gender, grade and stage.

Figure 5 Establishment and validation of the ICDRG risk models in HCC. (A) Forest plot based on univariate Cox analysis in the TCGA cohort. (B) Tenfold cross-validation for tuning parameter selection in the Least absolute shrinkage and selection operator (LASSO) model. Eleven ICDRGs expression were selected by the LASSO Cox models. (C) LASSO coefficients of the ICDRGs. (D) Nomogram for predicting 1-, 3-, and 5-year OS in patients with HCC. Time-dependent receiver operating characteristic (ROC) curve analysis in the TCGA (E) and ICGC (F) cohorts. ROC curves of clinical parameters in the TCGA (G) and ICGC (H) cohorts. (I) The calibration curves for 1-, 3-, and 5-year overall survival. Clinical impact curves for predicting OS in HCC patients in the TCGA (J) and ICGC (K) cohorts. Univariate (L) and multivariate (M) Cox regression analysis of clinicopathological features.

The Expression Levels of Prognosis-Related ICD Genes in HCC

To further elucidate the relationship between prognosis genes and HCC, we conducted qPCR analysis on a human normal liver cell line (LO2) and two HCC cell lines (HepG2 and Huh7) (Figure 6A). The results demonstrated that the expression of ATG5, BAX, CASP8, HSP90AA1, LY96, NT5E, P2RY2 and TLR2 was significantly upregulated in the liver cancer cells. These findings were consistent with results obtained from the TCGA cohort (Figure 6B). HPA database was performed to explore the protein expression of eight ICDRGs in the risk model. The protein expressions of LY96, NT5E, and TLR2 in HCC tissue were not included in the HPA database. We acquired the IHC results of the corresponding normal tissues and tumor tissues in the HPA database (Figure 6C). The results showed that ATG5, BAX, and CASP8 were highly expressed in HCC compared with normal liver tissue. However, there was no significant difference in the expression of HSP90AA1 and P2RY2 between normal liver tissues and HCC tissues.

Figure 6 The expression levels of prognosis-related ICD genes in HCC. (A) The expression levels of 8 prognosis genes in human normal liver cell line (LO2) and two HCC cell lines (HepG2 and Huh7) were examined by qRT-PCR. (B) The expression levels of 8 prognosis genes in the TCGA dataset. (C)The expression level of ICDRGs determined by immunohistochemistry in cancer tissues and normal tissues obtained from HPA datasets. (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

The Risk Score Was Related to HCC Prognosis and Pathological State

We classified all patients into high-risk and low-risk groups based on the median value of the risk score. As the risk score increased, patients had poorer OS and higher mortality in both the TCGA and ICGC cohorts (P<0.05, Figure 7A). The risk score and vital status among HCC patients were depicted by scatter plots and risk curves (Figure 7B). We found that the higher risk scores were significantly associated with ICD pathological grades and tumor stages in the TCGA cohort (Figure 7C–E). The relationship between survival status, subtypes and risk scores were further examined, and patients who died or had poor prognosis had higher risk scores (Figure 7F–G). Next, we confirmed the correlations between the risk scores and the expression of immune cells. As shown in Figure 7H-J, the risk scores were positively associated with Macrophages M0, but negatively associated with NK activated cells and CD8⁺T cells. The results showed that there were considerable differences in the ICDRG signatures among different risk groups (Figure 7K and L).

Figure 7 The risk score was related to HCC prognosis and pathological state. (A)K-M survival analysis of the ICDRGs risk model in the TCGA (left) and ICGC (right) cohorts for HCC patients. (B) Risk triple plots, including risk dispersion plots, survival time scatter plots, and heatmaps of model gene expression in the TCGA (left) and ICGC (right) cohorts. Boxplots of risk scores in HCC patients with different stages (C), T stages (D), pathological grades (E), status of survival (F), and clusters (G). (H-J) Correlation analysis between risk scores and immune cells. (K) Comparison of the differences in subtypes between different risk groups. (L) The relationship between ICDRGs and clinical indicators and subtypes of HCC.

Differences in Biological Functions Between Different Risk Groups of HCC

In order to identify the molecular mechanisms regulating prognosis, we further identified the key 2170 DEGs in high- and low-risk groups (Figure S2) and then performed GO and KEGG pathway enrichment analysis on the above DEGs. Based on the results of the GO analysis, the DEGs were predominantly enriched in pathways related to extracellular matrix organization, extracellular structure organization, negative regulation of hydrolase activity and wound healing (Figure 8A–C). Based on the results of the KEGG analysis, the DEGs were predominantly enriched in pathways related to tyrosine metabolism, glycolysis/gluconeogenesis, retinol metabolism, drug metabolism (cytochrome P450), and proteoglycans in cancer (Figure 8D–F). GSEA was further performed to complement and validate the functional annotation of KEGG and GO. According to GO enrichment analysis, the most enriched biological processes in the high-risk group were strongly connected to biological adhesion, external encapsulating structure organization, collagen containing extracellular matrix, external encapsulating structure, and extracellular matrix structural constituents (Figure 8G). Lipid metabolic processes, monocarboxylic acid metabolic processes, organic acid metabolic processes, small molecule biosynthetic processes, and small molecule metabolic processes were enriched in the low-risk groups (Figure 8H). According to KEGG enrichment analysis, the most enriched biological processes in the high-risk group were strongly connected to antigen processing and presentation, cell cycle, leishmania infection, oocyte meiosis and progesterone-mediated oocyte maturation (Figure 8I). Complement and coagulation cascades, drug metabolism cytochrome p450, metabolism of xenobiotics by cytochrome p450, retinol metabolism and steroid hormone biosynthesis were enriched in the low-risk groups (Figure 8J).

Figure 8 Differences in biological functions between different risk groups. The bar plot (A), bubble plot (B), and circular plot (C) of the GO pathways enrichment. The bubble plot (D), bar plot (E), and circular plot (F) of KEGG pathways enrichment. The gene set enrichment analysis (GSEA) for GO (G and H) and KEGG (I and J) for high-risk and low-risk groups of HCC patients.

Tumor Mutation Burden and Microenvironments Landscape in Different Risk Groups of HCC

The relationship between the prognostic characteristic, clusters, and clinical features was clearly shown in the Sankey diagram (Figure 9A). Most patients in C1 and C2 were in the low-risk group and had better prognostic outcomes, while the opposite results were observed in C3 and C4. Additionally, the majority of HLA genes and immunological checkpoints were elevated in the high-risk group (Figure 9B and C). On the contrary, the opposite trend was observed in the low-risk group. Patients in the high-risk group had significantly higher percentages of activated CD4⁺ T cells, activated dendritic cells, immature B cells, immature dendritic cells, MDSC, mast cells, natural killer T cells, neutrophils, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, type 2 T helper cells, and effector memory CD8⁺ T cells (Figure 9D and F). Figure 9E showed the mutation differences of the top 20 genes between different subtype groups in HCC. The Mutation of TP53 in the high-risk group had a higher level of mutation than in the low-risk group.

Figure 9 Tumor mutation burden and microenvironment landscape in different risk groups. (A) The Sankey diagram revealed the connection between cluster, risk score, and survival status. The boxplot illustrating the difference in human leukocyte antigen (HLA) (B) and immune checkpoints (C) and between the high-risk and low-risk groups. The boxplot (D) and heatmap (F) of immune infiltration cells between four subtypes of HCC. (E) Waterfall plots of gene mutations in low-risk (left) and high-risk (right) groups. (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

Chemotherapy Drug Sensitivity in Different Risk Groups of HCC Patients

We evaluated the IC50 of different chemotherapeutic drugs among different ICDRG risk groups in order to evaluate the potential use of ICDRGs in the individualized treatment of HCC. The results of the drug sensitivity showed that the low-risk group had lower IC50 of the following drugs than the high-risk group: topotecan, sorafenib, linsitinib, JQ1, and gemcitabine, which suggested that patients in the low-risk group displayed higher sensitivity from the chemotherapies mentioned above. While 5-fluorouracil, Paclitaxel, Fludarabine, Dasatinib, Cediranib, Vinorelbine, LCL161, Lapatinib, IWP-2, and Gefitinib had lower IC50 values in the high-risk group, suggested that patients in the high-risk group showed greater sensitivity to the chemotherapies mentioned above (Figure 10).

Figure 10 Chemotherapeutic response analysis. Boxplots comparing differences in half-inhibitory concentration (IC50) values between high- and low-risk score groups.

Discussion

The onset of HCC was insidious. The early symptoms were atypical or difficult to diagnose, and most were already in advanced stages when diagnosed.³⁰ The lack of safe and effective treatment was responsible for the progression of HCC and the increase in mortality.⁴ In this study, an ICD-related risk model consisting of eight genes was established to predict the prognosis of HCC patients. The results also revealed the diversity of the tumor immune microenvironments in the subtypes and risk groups, which also predicted the effectiveness of applied immunotherapy and chemotherapy, providing evidence for individualized treatment.

We determined 55 ICDRGs and analyzed their interactions, 35 of which were differentially expressed between HCC tissues and normal tissues. In the TME, ICD plays a major role in stimulating the dysfunctional antitumour immune system.³¹ ICD represents a particular form of cancer cell death accompanied by production of DAMPs that can be recognized by pattern recognition receptors, which enhances infiltration of effector T cells and potentiates anti-tumor immune responses.³² The abnormal functioning of ICD inevitably leads to an imbalance in the anti-tumor effect, consequently contributing to the occurrence and progression of HCC. The most direct cause of this imbalance is the differential expression of ICD-related genes. Therefore, we focused on the these 35 genes when we built the risk model next.

There was a significant association between ICDRGs and prognosis. By subtype analysis, we found that patients with low expression levels of ICDRGs had a better prognosis. In contrast, patients with high expression levels of ICDRGs had a poor prognosis. In addition, the patients with higher expression of ICDRGs were accompanied by tumor heterogeneity and extensive genomic alterations, primarily in the form of high-frequency TP53 mutations, higher stromal and immune scores, and lower tumor purity. Aran et al³³ suggested that tumor samples with lower tumor purity contain more immune cells, allowing for an increased inflammatory response and high mutation load in tumor cells. This corresponded to the results we observed. Immune infiltration analysis showed that the subtypes with the better prognosis had more innate immune cell expression. Li et al³⁴ suggested that triggering ICD in HCC cells facilitated the ability of dendritic cells and macrophages to recognize and phagocytose tumor cells effectively, which ultimately resulted in the activation of an anti-tumor immune response. It might help to explain our findings. We also found that the amount of immune cells with phagocytic and cytotoxic functions (such as macrophages and NK cells) was reduced in the subtype with poor prognosis. These results pointed to the possibility that ICDRG might lead to different prognostic outcomes by modulating immune cells.

In order to provide a better application of ICDRGs to the clinical treatment of HCC, we developed a new prognostic model that includes eight genes associated with ICD (ATG5, BAX, CASP8, HSP90AA1, LY96, NT5E, P2RY2 and TLR2) and created a nomogram containing clinical indicators. The model demonstrated good discrimination, calibration and prognosis-predictive capacity. In the ICGC cohort, we successfully carried out an external validation of the model.

Previous studies demonstrated that these ICDRGs were strongly associated with the development and prognosis of cancer. It was demonstrated that suppressing ATG5 induced autophagy inhibition and overcame drug resistance, thereby enhancing the growth inhibition of HCC by chemotherapeutic drugs.^35–37 Bax regulated HCC cell proliferation and radiosensitivity, and BAX redistribution induced apoptosis resistance and selective stress sensitivity.^38,39 Studies showed that Casp-8 upregulation significantly reduced HCC programmed cell necrosis and increased resistance to apoptosis, thereby promoting the malignant progression of HCC.^40,41 P2Y2R stimulated the DNA damage response and hepatocyte proliferation, therefore promoting hepatocarcinogenesis.⁴² Previous research has shown that TLR2-dependent autophagy induced by hepatoma-derived factors promotes M2 macrophage differentiation, reinforcing tumor malignant behaviors.⁴³ These findings backed up our immunohistochemical findings that the majority of these genes were elevated in HCC tissue and were associated with a poor prognosis.

We used the ICDRG risk model to calculate the risk score for each HCC patient. The results showed that the risk score was an independent risk factor for a poor prognosis. The risk scores were positively correlated with the severity of the pathology (tumor stages and pathologic grades) but negatively correlated with immune cells with killing function, such as NK cells and CD8⁺ T cells. High- and low-risk groups were distinguished by the risk scores. The K-M survival analysis and risk curve analysis showed a poor prognosis for the high-risk group. Multiple bioinformatics analysis revealed significant differences between high- and low-risk groups in terms of genetic landscape, immunological status, and pharmaceutical sensitivity. Immune infiltration analysis revealed that patients in the high-risk group recruited a variety of immune-suppressive cells, including Tregs, Th17 cells, and MDSCs. The above evidences suggested that the suppressive immune microenvironments and tumor immune escape were more common in high-risk HCC patients, resulting in a poor prognosis.

Human leukocyte antigen-I (HLA-I) played an important role in antitumor immune response and tumor progression. Patients with less HLA diversity and fewer tumor mutations respond poorly to immune checkpoint inhibitor treatment.⁴⁴ Interestingly, almost all the suppressive immune checkpoints and HLA genes were significantly upregulated in the high-risk group. The results described above supported that immune checkpoint inhibitor therapy might be beneficial for high-risk patients and also confirmed the applicability of the ICDRG risk model for predicting immunotherapy. Notably, due to drug resistance and tumor heterogeneity, HCC patients in different risk groups responded differently to various chemotherapeutic drugs. Our results demonstrated that the ICDRG risk model was able to predict the sensitivity of chemotherapeutic drugs for HCC patients based on the different inhibitory concentrations (IC50) observed in different risk groups.

Our study still had some limitations. First, since this study was a retrospective investigation, prospective studies with real-world analysis are necessary to validate the application of this strategy. Second, cell-based validation is necessary for the tumor immune cell landscape within the tumors of HCC patients that was identified using gene expression data. Third, the absence of animal experiments restricts our ability to directly assess the translational applicability of our findings in vivo.

Conclusion

This study identified ICDRGs in HCC patients, developed, and validated an ICDRGs risk model to predict OS in HCC patients, demonstrating good predictive performance. We also investigated the differences in immunotherapy response and chemotherapeutic drug sensitivity between risk groups, which were associated with genomic heterogeneity and immune-suppressive status. Our study provided a new perspective for assessing the prognosis of HCC patients and new individualized treatment strategies.