Abstract Accumulating studies have highlighted the biological significance of immunogenic cell death (ICD) in cancer immunity. However, the influence of ICD on tumor microenvironment (TME) formation and immune response in Hepatocellular carcinoma (HCC) remains largely unexplored. In this study, we systematically analyzed the mRNA profiles of ICD-related genes in 1847 HCC patients and identified three molecular subtypes with significantly different immune features and prognostic stratification. A reliable risk model named ICD score was constructed via machine learning algorithms to assess the immunological status, therapeutic responses, and clinical outcomes of individual HCC patients. High ICD score indicated an immune-excluded TME phenotype, with lower anticancer immunity and shorter survival time. In contrast, low ICD score corresponded to abundant immune cell infiltration, high sensitivity to immunotherapy and a positive prognosis, indicating an “immune-hot” phenotype. Pan-cancer analysis further validated a negative association between ICD score and the immune cell infiltration levels. In conclusion, our findings revealed that the ICD score could serve as a robust prognostic biomarker to predict the benefits of immunotherapy and optimize the clinical decision-making of HCC patients. Keywords: Immunogenic cell death, Tumor microenvironment, Prognosis, Precision treatment, Hepatocellular carcinoma 1. Introduction Primary liver cancer is the sixth most prevalent malignancy and the second leading cause of cancer related deaths worldwide [[39]1,[40]2]. Hepatocellular carcinoma (HCC) accounts for approximately 90 % of liver cancer cases and is expected to affect over one million people by 2025 [[41]3]. Although surgical resection, ablation, and liver transplantation may offer favorable outcomes, most HCC patients are diagnosed at advanced stages that are usually ineligible for curative treatments [[42]4]. Systemic chemotherapy remains a valuable option for these patients [[43]5], but limited survival benefit, drug resistance and adverse side effects greatly hinder the clinical benefits [[44]6]. Immune checkpoint inhibitors (ICIs) have revolutionized the therapeutic outcomes of advanced cancers. Unfortunately, only a small fraction of HCC patients have gained considerable benefits [[45][7], [46][8], [47][9]]. Therefore, identifying reliable prognostic biomarkers is essential for optimizing drug therapies and enhancing the survival of HCC patients in the era of individualized treatment. As a typical inflammation-associated tumorigenesis, the immunosuppressive microenvironment of HCC is closely correlated with immune tolerance and evasion [[48]10]. Immunogenic cell death (ICD) is a distinct form of immunoregulatory cell death (RCD), which can activate the tumor microenvironment (TME) by inducing ICD in tumor cells [[49]11]. It is characterized by the release of danger-associated molecular patterns (DAMPs) and tumor-associated antigens (TAAs) [[50]12]. Surface-exposed CRT binds to the CD91 receptor, serving as a key “eat me” signal that promotes phagocytosis by antigen presenting cells (APCs) [[51]13]. Extracellular ATP released from dying cells operates as a “find me” signal to recruit phagocytes such as dendritic cells (DCs), monocytes, and macrophages into the tumor site [[52]14]. The nuclear chromatin binding protein HMGB1 interacts with Toll-like receptor 4 (TLR-4), promoting DCs maturation and stimulating proinflammatory cytokine release [[53]15]. Intracellular HSPs including HSP70, are highly conserved proteins that attract phagocytes and facilitate the activation of natural killer (NK) cells [[54]16]. DAMPs enhance the recruitment and maturation of APCs, facilitate TAA presentation to effector T cells, activate cytotoxic T lymphocytes (CTLs), and promote the secretion of immunostimulatory cytokines, ultimately stimulating tumor-specific immune responses and provoking long-term anticancer immunity [[55]17]. These studies highlight the fundamental roles of ICD in the TME and cancer immunotherapy. A deeper understanding of the correlation between ICD and TME cells may point us to new biomarkers that can help optimize precision treatment and the development of immunotherapy for patients with HCC. In this research, we systematically investigated the expression profile of ICD-related genes and identified three molecular subtypes in 1847 HCC patients from 9 independent public datasets. The differences between three ICD patterns in functional annotations, prognostic status and immunological features were compared. Using machine learning algorithms, we established a risk stratification signature to predict prognosis, immune characteristics, and immunotherapy efficacy in HCC patients. Finally, a comprehensive pan-cancer analysis was conducted to explore the correlation between ICD score, clinical outcomes, and immune cell infiltration across various cancers. 2. Materials and methods 2.1. Data collection and processing The overall workflow of our research is illustrated in [56]Fig. S1. Pan-cancer transcriptome profiling, clinicopathological information, copy number variation (CNV), somatic mutation, and methylation data were sourced from the UCSC Xena Database. We systematically searched the clinical features and gene expression data of HCC samples in The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and International Cancer Genome Consortium (ICGC). Patients who fulfilled the following selection criteria were obtained: (a) histologically diagnosed liver cancer; (b) detailed RNA sequence data; (c) cohorts >50 cases; (d) delete technical replications if necessary. Ultimately, a sum of 1847 HCC patients from 9 independent public datasets were included for further analysis. Among them, four datasets (TCGA-LIHC, ICGC-LIHC, [57]GSE76427 and [58]GSE14520) with complete clinical information and overall survival (OS) data were combined into one meta-cohort. The fragments per kilobase per million (FPKM) were converted to transcripts per million (TPM). CHCC-HBV (Chinese HCC patients with hepatitis B virus infection) cohort data was obtained from NODE and the supplementary materials of Gao et al. [[59]18]. The Clinical Proteomic Tumor Analysis Consortium (CPTAC) website was applied to extract the gene-level proteomics data of CHCC. Four datasets ([60]GSE6764, [61]GSE69164, [62]GSE77509 and [63]GSE109211) were obtained to assess the tumor progression, portal vein tumor thrombus (PVTT) presence and sorafenib sensitivity. In addition, we enrolled [64]GSE35640, [65]GSE91061 and Nathanson's cohort [[66]19] to assess the ability of ICD score in predicting immunotherapy response. 2.2. Consensus clustering Unsupervised clustering analysis was conducted to determine ICD phenotypes in HCC. The “ConsensusClusterPlus” package was used to identify the optimum cluster numbers, with 1000 repetitions performed to ensure classification stability [[67]20]. 2.3. Estimation of immune cell infiltration The immune infiltration was evaluated by eight independent algorithms including ESTIMATE, Timer, EPIC, MCP-counter, xCell, Quantiseq, CIBERSORT and CIBERSORT-abs. We applied the ssGSEA algorithm to evaluate the enrichment scores for the seven critical steps of cancer immunity cycles [[68]21]. 2.4. Functional enrichment The correlations between ten cancer related pathways and ICD-related molecules were explored. The pathway scores were determined by summing the relative protein levels of all positive regulatory components and subtracting those of the negative regulatory components [[69]22]. The gene sets in “h.all.v7.4.symbols.gmt” obtained from the MSigDB database were applied for Gene Set Variation Analysis (GSVA) analysis [[70]23]. 2.5. Single-cell RNA-seq (scRNA-seq) analysis Two scRNAseq datasets ([71]GSE149614 [[72]24] and [73]GSE125449 [[74]25]) were obtained from GEO database and further reanalyzed in this study. [75]GSE149614 recruited 10 HCC patients with scRNA-seq data of four relevant tissue types. Furthermore, the single-cell transcriptome data of 9 HCC tumors were downloaded from [76]GSE125449 dataset. 2.6. Somatic mutation and copy number alteration analysis The R package “TCGAbiolinks” was used to download somatic mutation and CNV data of TCGA-LIHC from Genomic Data Commons (GDC) [[77]26]. The gene copy number amplification/deletion data and GISTIC score were calculated using GISTIC2.0 in GenePattern. The individual fraction genome altered (FGA) value was analyzed by methods from previous research [[78]27]. 2.7. Drug sensitivity estimation The Broad Institute-Cancer Cell Line Encyclopedia (CCLE) project was applied to obtain the transcriptome profiling data of cancer cell lines (CCLs) [[79]28]. The CERES scores of genome-wide CRISPR knockout screens were retrieved from DepMap (Cancer Dependency Map). CERES algorithm was used to assess the dependency of candidate genes, and the gene with lower score indicated greater influence on proliferation and survival [[80]29]. The R package “pRRophetic” utilized drug sensitivity data from the Cancer Treatment Response Portal (CTRP) and PRISM Reuse datasets to predict chemotherapy sensitivity for each HCC patient [[81]30]. We retrieved the half-maximal inhibitory concentration (IC50) of antitumor drugs from the Genomics of Drug Sensitivity in Cancer (GDSC) database, calculated the IC50 value for each sample, and applied 10-fold cross-validation to evaluate the prediction accuracy [[82]31]. Lower IC50 (or AUC) value indicated increased sensitivity to specific compound. Using the connectivity map (CMap) data, the XSum method was performed to forecast the possible small-molecule drugs for HCC [[83]32,[84]33]. 2.8. Construction of the risk score model The random forest (RF) algorithm was utilized to assess the importance of 27 ICD-related molecules in prognosis, and genes with significance scores >0 were screened. The classification of gene combinations was identified by Gaussian mixture model (GMM)-based hierarchical clustering method. Genes in the clusters with highest AUC value were subjected to the least absolute shrinkage and selection operator (LASSO) analysis. Then, the ICD score was established using the Cox coefficients (β) of ten genes identified by LASSO regularization as follows: [MATH: Riskscore=∑i=1n Geneexpr i×βi :MATH] Patients were dichotomized into high- and low-risk groups using optimal cut off value identified by the R package “survminer”. 2.9. Statistical analyses All statistical analyses were conducted using R 4.2.2 software. Pearson or Spearman correlation analysis was applied to estimate the association between two variables. The differences of continuous variables between two or more groups were compared by Wilcoxon signed-rank test and Kruskal-Wallis test. Chi-square test was used to test the categorical features of different subtypes. All P-values were two-sided and P < 0.05 was considered statistically significant. 3. Results 3.1. Multi-omic alterations of 27 ICD-related molecules We conducted a comprehensive analysis of the genomic variation, gene expression and DNA methylation profiles of 27 ICD-related molecules in the TCGA pan-cancer database. The waterfall plot was used to visualize the top 10 mutated genes. Of the 3146 samples, 2350 exhibited mutations of ICD genes with a frequency of 74.4 % and missense mutation was the most prevalent type ([85]Fig. 1A). PIK3CA exhibited the highest mutation frequency (44 %) among the mutated genes, followed by NLRP3 and TLR4, while HMGB1 and ATG5 showed relatively low mutation rate ([86]Fig. 1A and B). Additionally, we found that the THYM and PCPG cohorts showed relatively low mutation frequency, while the overall mutation rate was relatively high in BRCA and UCEC cohorts. Comparing gene methylation differences between tumor and normal samples revealed complex patterns of ICD-related molecules across 14 cancer types ([87]Fig. 1C). The gene expression levels and DNA methylation were negatively correlated in most cancers ([88]Fig. 1D). The CNV percentage analysis indicated that NLRP3 and PIK3CA exhibited the highest frequency of CNV, and heterozygous amplification and deletion were the main CNV types ([89]Figs. S2A–C). We then examined the impact of CNA on the gene expression and found varying associations between mRNA expression of ICD-related genes and CNA across different cancer types ([90]Fig. 1E). Consequently, our findings revealed that genome variance can significantly regulate the expression patterns of ICD-related genes. Fig. 1. [91]Fig. 1 [92]Open in a new tab Expression variation of ICD-related molecules. (A) The waterfall diagram showing the 10 most frequency mutated genes across pan-cancer in the TCGA cohort. (B) The mutation frequency of ICD-related molecules in different cancer types. Numbers represent the number of samples with corresponding mutated genes. (C) The bubble diagram showing the differences in methylation expression of ICD-related molecules between cancer samples and normal samples. (D) Correlation between methylation and mRNA expression level. (E) Correlation between CNV and mRNA expression of the ICD-related genes. (F) The histogram (upper panel) displaying the number of significantly differentially expressed genes, and the bubble chart shows the fold change (FC) and FDR of ICD-related genes in each cancer. (G) Summary of the association between expression of ICD-related molecules and patient survival in each cancer type. (H) The location of CNV alteration on chromosomes in the TCGA-HCC cohort. (I) The frequency of CNV gain (red), loss (blue) and non_CNV (green) of ICD-related molecules. We analyzed 20 cancer types to compare ICD expression profiles in primary tumors and adjacent normal tissue. ICD-related molecules were extensively dysregulated in various cancer types compared to their normal samples ([93]Fig. 1F). Among the significantly differentially expressed genes, the majority were downregulated in LIHC, LUAD, and LUSC, but upregulated in GBM, ESCA, KIRC, HNSC, and KIRP, indicating that ICD might play different roles in different cancers. To further explore the prognostic significance of ICD-related genes, we conducted a pan-cancer survival analysis and found that these genes exhibited strong correlation with clinical outcomes in several types of cancer ([94]Fig. 1G). To investigate the genomic features of 27 ICD-related molecules in HCC cohort, the rates of somatic mutations and copy number variations were summarized. The DNA mutation of ICD-related genes occurred in about 17.31 % of patients. Among them, PIK3CA experienced the highest mutation frequency of up to 3 %, with missense mutations being the predominant type ([95]Fig. S2D). The location of CNV alteration of ICD-related genes on chromosomes was depicted in [96]Fig. 1H. In terms of somatic copy number alterations, NLRP3 showed highest frequency of CNV gain, followed by PIK3CA ([97]Fig. 1I). 3.2. Identification of ICD subtypes in HCC In order to identify the regulatory mechanism of ICD-related molecules in HCC, we included a total of 1847 samples obtained from nine cohorts. The baseline of each cohort was depicted before ([98]Fig. 2A) and after ([99]Fig. 2B) the batch effect which conducted by “ComBat” algorithm. It can be seen that we had effectively corrected the batch effects between different datasets. Spearman correlation analysis revealed significant associations between the 27 ICD-related genes and tumor-infiltrating immune cells ([100]Fig. S3A). The prognostic values of ICD-related genes in patients with HCC were evaluated by univariate Cox regression model. We found that ATG5, BAX, CALR, EIF2AK3, HMGB1, HSP90AA1 and PDIA3 were related to a poor prognosis, whereas CD4 and CD8A appeared to be associated with a superior prognosis ([101]Figure S 3B). We subsequently utilized a network diagram to comprehensively depict the landscape of 27 related molecules interactions and their prognostic significance ([102]Fig. 2C). The expression levels of ICD-related genes exhibited significantly positive correlations both within the same category and across different modification patterns, indicating that the imbalance of these molecules was strongly linked to HCC occurrence and progression. Fig. 2. [103]Fig. 2 [104]Open in a new tab Unsupervised learning to determine three ICD phenotypes. (A–B) Principal component analysis (PCA) showing the gene expression distribution of ICD-related molecules in nine HCC cohorts before (A) and after (B) batch effect correction. (C) The interaction between ICD-related molecules based on the meta-cohort. The circle size represents the survival impact of each gene and scaled by P-value. The thickness of the line indicates the correlation strength assessed by Spearman correlation analysis. (D) Histogram showing the distribution of ICD subtypes in nine HCC cohorts. (E) Kaplan-Meier curves of OS for HCC patients in the meta-cohort with three ICD clusters. (F) Differences in the expression of ICD-related molecules among three subtypes. Statistical significance is shown by asterisks (∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001; ∗∗∗∗P < 0.0001). (G) Differences in the enrichment scores of cancer immunity cycles among three subtypes. Fig. 3. [105]Fig. 3 [106]Open in a new tab Regulation of immunomodulators (IMs) and analysis of biological functions in distinct ICD clusters. (A) From left to right: mRNA expression (median normalized expression levels); expression vs. methylation (the association between gene expression and DNA-methylation beta-value); amplification frequency (the difference between the fraction of samples in which an IM is amplified in a specific subtype and the amplification fraction in all samples); and the deletion frequency (as amplifications) for IM genes by ICD clusters. (B) The heatmap displaying the correlation between the mRNA expressions of ICD-related genes in cancer-related signaling pathways. Only a function (inhibit or activate) in at least 5 cancer types was shown. (C) The bubble chart showing the differences in enrichment scores of carcinogen pathways among three clusters. (D) Enrichment of cancer hallmark signaling pathways in different subtypes. (E) Score variations in biological pathways among three ICD subtypes. Three distinct modification patterns were identified based on the transcriptomic profiles of 27 ICD-related genes using the unsupervised clustering method ([107]Fig. 2D, [108]Fig. S3C), named Cluster A (748 cases), Cluster B (296 cases), and Cluster C (803 cases). Survival analysis of the meta-cohort with complete clinical information indicated that patients in different clusters exhibited significantly different survival rates, and Cluster C had the worst prognosis (log-rank, P = 0.03, [109]Fig. 2E). The mRNA expression levels of 27-related genes varied significantly among the three clusters ([110]Fig. 2F). Specifically, Cluster B demonstrated the highest expression levels of ICD-related molecules, followed by Cluster A and Cluster C. Six algorithms were applied to assess the tumor immune infiltration in HCC. We found that Cluster B had significantly higher quality immune cells than the other subgroups ([111]Fig. S3D). Similarly, gene set enrichment analysis revealed significantly enhanced cancer immunity cycle activities in cluster B ([112]Fig. 2G). 3.3. Molecular characteristics and biological functions in different ICD subtypes Given the importance of immunomodulators (IMs) in cancer immunotherapy, it is essential to gain a further understanding of their expression in different ICD subtypes [[113]34]. Our findings indicated that the expression levels of IM genes varied across ICD subgroups, with significantly higher expression observed in cluster B compared to the other clusters ([114]Fig. 3A). The correlation between DNA methylation and IM genes also differed among different ICD groups. Notably, several genes including CD28, VEGFB, CD27, CD40, ADORA2A and ARG1 demonstrated positive associations with DNA methylation, indicating epigenetic silencing. CNVs affected many IMs and varied across three ICD subtypes. In particular, Cluster A and B were found to exhibit frequent amplification or deletion of most IM genes. We then investigated the characteristics of several carcinogen and immune-related signaling pathways. As depicted in [115]Fig. 3B, ICD-related genes were negatively associated with DNA damage response, cell cycle, and hormone AR activation, while consistently promoting EMT, apoptosis, and hormone ER signaling in pan-cancer. Previously published signatures were applied to further examine the differences in enrichment scores of carcinogen pathways across three clusters [[116]35]. We found that Cluster A showed higher enrichment scores for DNA repair, while Cluster B was significantly enriched in the PI3K-AKT-mTOR and Wnt/β-Catenin pathways ([117]Fig. 3C). Cluster C showed significant association with MYC and TGF-β activation. Through GSVA analysis with cancer hallmark gene sets, we observed that cluster B was significantly associated with inflammatory signaling pathways, whereas cluster C was linked to metabolism-related signatures ([118]Fig. 3D). Additionally, we conducted a comparison of the average pathway scores between different subtypes and found that cluster A was characterized by high levels of Homologous recombination (HR), Mismatch repair (MMR), Base excision repair (BER), and Nucleotide excision repair (NER), while cluster B exhibited high levels of CD8 T-effector, TMEscore, and RAS ([119]Fig. 3E). To decipher the mutational features in three ICD clusters, we analyzed the genomic variation data derived from TCGA cohort. Cluster A exhibited more mutations of HMCN1, CACNA1E and ABCA13, while Cluster C was enriched for mutations in OBSCN, AXIN1, and ALB. Cluster B exhibited a lower mutation frequency than Clusters A and C ([120]Fig. S4A). No statistically significant difference was found in tumor mutation burden among the three subtypes ([121]Fig. S4B); however, the median number of TMB for Cluster A was significantly higher compared to Cluster C. Furthermore, Cluster B showed lower copy number alterations (CNAs), indicating a better chromosomal stability ([122]Figs. S4C and D). 3.4. Establishment of the ICD-related prognostic signature A prognostic model named ICD score was developed to evaluate the ICD modification pattern in HCC patients. The importance of 27 genes was determined by the random forest algorithm, and 13 genes with importance scores greater than 0 were finally selected ([123]Fig. 4A and B). As depicted in [124]Fig. 4C, the remaining 13 candidates had 8191 combinations. We performed GMM and ROC analysis to cluster these gene sets, and found that combinations of Cluster 1 which contained 10 genes showed the highest AUC value. Afterward, we utilized the LASSO Cox regression analysis to create our risk score model ([125]Figs. S5A and B), and the ICD score was then calculated for each HCC patient. Using the optimal cut-off value from the “survminer” package, HCC patients were dichotomized into high- and low-risk groups. Fig. 4. [126]Fig. 4 [127]Open in a new tab Construction and validation of the ICD score. (A) Correlation between the number of random forest trees and error rates. (B) The importance of selected features based on the RF model. (C) Different GMMs based on the expression values of 13 candidates using clustering analysis. (D) Kaplan-Meier survival curves for low- and high-risk groups in the meta-cohort. (E) The distribution of ICD score, survival status, and the ten-gene expression panel. (F) Univariable and multivariable Cox regression analysis of ICD score and prognosis in the meta-cohort. (G–I) Kaplan-Meier survival curves (upper panel) for low- and high-risk groups, and the correlation between ICD score and clinicalpathological characteristics in TCGA (G), ICGC (H) and [128]GSE14520 (I) cohorts. We found that patients with high ICD scores showed significantly dismal overall survival compared to those with low ICD scores in the meta-cohort ([129]Fig. 4D). The distribution of the ICD scores, overall survival of HCC patients, and expression profiles of selected genes were illustrated in [130]Fig. 4E. We sought to identify whether ICD score could be an independent prognostic biomarker for HCC patients. The statistical significance of the model for OS was confirmed by Univariate and multivariate Cox regression analyses, demonstrating that ICD score was a reliable biomarker to predict clinical outcome ([131]Fig. 4F). Next, we validated the reliability of ICD score in five independent HCC cohorts including TCGA, ICGC, [132]GSE14520, [133]GSE76427 and CHCC-HBV cohorts. Consistent with the above findings, high-risk score was obviously associated with unfavorable outcomes ([134]Fig. 4G–I, Figure S5C, D). Additionally, our research revealed that ICD score was obviously associated with different clinicopathological features. The ability of ICD score to distinguish HCC progression was further examined in [135]GSE6764 cohort. Our findings demonstrated a gradual increase of ICD score in the process of Normal-very Early HCC-very Advanced HCC ([136]Fig. S5E). Furthermore, we included [137]GSE69164 and [138]GSE77509 cohorts to explore the correlation between ICD score and PVTT. The results revealed that the risk score of adjacent normal tissue was significantly lower than that of PVTT ([139]Figs. S5F and G). 3.5. Functional enrichment in high- and low-risk groups We calculated the enrichment scores of the hallmark pathway gene signatures to investigate the differences in pathway enrichment analysis between the two subtypes. Direct comparison of high-versus low-risk score group revealed significantly distinct biological functions ([140]Fig. 5A). The high-risk group showed significant enrichment in stromal and carcinogenic activation pathways, whereas the low-risk group was enriched in immune activation and metabolism-related biological process. Subsequent analyses suggested that the high-risk score group was prominently correlated with unfavorable molecular features such as MYC and TGF-β pathways ([141]Fig. 5B). To assess the stem cell-like properties in individual patients, we calculated the enrichment scores of stemness cell signatures ([142]Fig. S6A). The ICD score was positively correlated with stemness, indicating that high-risk score group tended to exhibit reduced progenitor-like characteristics and proliferative capacity. Fig. 5. [143]Fig. 5 [144]Open in a new tab Clinical significance and immune landscape of two ICD subtypes in the meta-cohort. (A) GSVA enrichment analysis showing the activation status of biological pathways in high- and low-risk group. (B) Differences in the enrichment scores of tumor-related signaling pathways between the two ICD groups. (C) Comparisons of immune score, stromal score, ESTIMATE score and tumor purity between two subtypes. (D) Correlation between ICD score and the abundances of immune cells. (E) Correlation between ICD score and the activation of cancer immunity cycles. 3.6. Immune landscape and genomic variations in high- and low-risk groups Given the importance role of tumor immune microenvironment in regulating cancer cell fate, we attempted to explore the landscape of immune cell infiltration in two different ICD subgroups. We found that ICD score was negatively linked to the Immune, Stromal, and ESTIMATE scores, but positively related to tumor purity ([145]Fig. 5C). ssGSEA analyses revealed statistically significant differences in immune cell infiltration and immune function between the two subgroups ([146]Fig. S6B). Five algorithms were also applied to examine the robustness and stability of the above findings. The ICD score was inversely associated with the infiltration level of most immune cells, except for Neutrophils (MCP-counter), Mast cells (xCell), and Macrophages M1(xCell) ([147]Fig. 5D). Given these observations, we subsequently correlated the risk score with the critical steps of the cancer-immunity cycle. As predicted, high-risk group appeared to have lower activities of most anti-cancer immune responses, prompting a significant correlation between higher ICD score and immunosuppression in HCC patients ([148]Fig. 5E). We then analyzed two scRNA-seq datasets ([149]GSE149614 and [150]GSE125449) to highlight the pivotal role of ICD score in the TME of HCC. After quality control and removal of batch effects, we calculated the ICD score of four tissue types from [151]GSE149614 dataset. It was observed that normal tissue sample exhibited significantly lower ICD score compared to primary tumor, MLN and PVTT tissues (Kruskal−Wallis, P < 2.2e−16, [152]Fig. S7A). As for primary tumor tissues, qualified single cells were clustered and categorized into six cell types: myeloid cells, fibroblast cells, hepatocyte cells, endothelial cells, T cells and B cells ([153]Figs. S7B and C). We found that the ICD scores of hepatocyte, fibroblast, and endothelial cells were higher than that of T and B cells ([154]Figs. S7D and E). Furthermore, immune cells exhibited significantly lower ICD scores than malignant and stromal cells ([155]Fig. S6F). The expression patterns of ten selected ICD-related genes for all cell types were also visualized ([156]Figs. S8A and B). Of these genes, CALR, HMGB1, HSP90AA1, and PDIA3 were expressed in both immune and stromal cells, while ATG5, NLRP3 and TNF were detected at low level in all cell types. The scRNA-seq data of nine HCC tumors in [157]GSE125449 was analyzed to confirm the influence of ICD score on immune regulation. Consistently, we observed negative association between ICD score and the immune cell infiltration like CD4^+ T cell, CD8^+ T cell and B cell ([158]Figs. S9 and S10). Taken together, these findings suggested that our ICD score was strongly associated with the formation of TME complexity and diversity. Finally, the association between ICD score and the genomic variations was examined. No significant difference in somatic mutations was found between the high-risk and low-risk groups ([159]Fig. S11A). The burden of copy number amplifications and deletions were revealed to correlate positively with ICD score at the arm-level ([160]Fig. S11B). We did not observe any significant difference at the focal level between the two subtypes. Interestingly, the levels of FGA, FGG, and FGL values were statistically higher in high-risk group and advanced HCC grade ([161]Fig. S11C). In addition, we intuitively displayed the topography of CNV in two subtypes. The high-risk group appeared to exhibit increased GISTIC score and copy number gain/loss percentage ([162]Fig. S11D). 3.7. Prediction of potential effective drugs Given that proteins highly correlated with ICD score might exhibit therapeutic significance for high-risk group, we attempted to identify druggable therapeutic targets based on the data obtained from CHCC-HBV cohort. We first calculated the correlation coefficient between ICD score and the expression level of druggable proteins. 216 protein targets with correlation coefficient greater than 0.20 were screened (P < 0.05, [163]Fig. 6A). We identified 29 poor prognosis-dependent targets by analyzing the correlation between CERES and ICD score in HCC cell lines (Spearman's r < −0.3 and P < 0.05, [164]Fig. 6B). The above two analyses ultimately yielded six genes including SOAT1, RPS6KA3, PTPN1, ACTR3, CHD1, and PDE6D. The CERES scores of RPS6KA3, PTPN1, ACTR3, CHD1, and PDE6D were generally less than zero in most hepatoma cell lines, suggesting their potential utility as therapeutic targets in HCC patients with high ICD scores. Fig. 6. [165]Fig. 6 [166]Open in a new tab Identification of potential therapeutic agents for HCC patients. (A) Volcano plot (left) and scatter plots (right) of correlation analysis and significance between CERES score of drug targets and ICD score. (B) Volcano plot (left) and scatter plots (right) of correlation analysis and significance between protein expression of drug targets and ICD score. (C–D) Drugs sensitive to chemotherapeutic agents in high-risk group based on the CTRP (C) and PRISM (D) datasets. (E) Top ranked six compounds with highest reversal potency were identified by the optimal method (XSum). (F–G) Drug sensitivity of Erlotinib (F) and Sorafenib (G) in two ICD subtypes. (H) Differences in the expression of ICD score between distinct clinical outcomes of Sorafenib. (I) Proportion of patients responding or not responding to Sorafenib in high- and low-risk groups. We adopted CTRP and PRISM databases to identify highly sensitive drug candidates for high-risk HCC patients. The association between ICD score and AUC value was assessed for each TCGA sample, and compounds with negative correlation coefficients (Spearman r < −0.23 for CTRP or −0.15 for PRISM) were further selected. We identified six CTRP- and PRISM-derived compounds with lower estimated AUC values that might be sensitive for patients with high ICD scores ([167]Fig. 6C and D). The similarity scores of all compounds obtained from CMap data were calculated using XSum method and a compromising parameter (topN = 500). Lower XSum scores indicated higher reversal potency and greater application potential. Six potential compounds with highest reversal potency that might have chemopreventive effects for HCC were identified, including X4.5.dianilinophthalimide, Nercaptopurine and Clofibrate ([168]Fig. 6E). Subsequently, we examined whether our ICD score could predict patients' response to receptor tyrosine kinase (RTK) inhibitors based on the GDSC database. It was observed that Erlotinib and Sorafenib exhibited significantly lower IC50 values in the low-risk group, indicating greater sensitivity for patients with low ICD scores ([169]Fig. 6F and G). The ICD score's predictive value in HCC patients undergoing Sorafenib treatment was further investigated in [170]GSE109211 cohort. We found that patients responding to Sorafenib generally had lower ICD scores, whereas non-responders were more prevalent in the high-risk subgroup ([171]Fig. 6H and I). To conclude, patients with low ICD scores were more likely to benefit from RTK inhibitors, particularly Sorafenib. 3.8. The role of ICD score in predicting immunotherapeutic benefits To explore the effect of ICD score on immune checkpoint blockade (ICB) response, we applied tumor immune dysfunction and exclusion (TIDE) to investigate the potential benefit from immunotherapy in different groups [[172]36,[173]37]. A high TIDE score was linked to immune escape, poor immunotherapy response, and unfavorable survival outcomes. We found that TIDE score was obviously elevated in the high-risk subtype, indicating that patients with high ICD scores were less likely to benefit from ICB therapy ([174]Fig. 7A). In particular, low-risk subtype had extensively higher levels of IPS z-score, MHC, and T cell dysfunction, but lower T cell exclusion compared to the high-risk subtype. No significant difference in suppressor cells was observed between the two subtypes ([175]Fig. 7A and B). Given the essential role of CTL in immune regulation and anti-cancer immunity, we calculated the CTL score by averaging the expression of five signature genes (CD8A, CD8B, GZMA, GZMB, and PRF1) to quantify the cytotoxicity of tumor-infiltrating T cells [[176]38]. Our analysis revealed that ICD score was negatively related to CTL score. The tumor inflammation signature (TIS) score is a quantitative index used to determine the effectiveness of cancer immunotherapy [[177]39]. Interestingly, we found that patients with high ICD scores showed a significantly lower TIS score than those with low ICD scores. These results indicated that low ICD score group appeared to have a better response to immunotherapy. Fig. 7. [178]Fig. 7 [179]Open in a new tab Prediction of response to immunotherapy. (A) Violin illustration showing comparisons of TIDE, Dysfunction, Exclusion and AZ values in high- and low-risk subtypes. (B) Violin illustration showing comparisons of MHC, CP, CTL and TIS values in high- and low-risk subtypes. (C, D) Bar and Violin graphs showing the treatment response to anti-CTLA4 therapy in high- and low-risk groups in Nathanson's cohort. (E) Kaplan-Meier survival curves for low- and high-risk groups in Nathanson's cohort. (F–I) Bar and Violin graphs showing the treatment response to immune checkpoint blockade in high- and low-risk groups in [180]GSE91061 and [181]GSE35640 database. (J) The submap algorithm for predicting the sensitivity to anti-PD1 and anti-CTLA4 treatment in patients with high and low ICD scores. Given the strong relationship between ICD score and immune response, we tried to elucidate whether our risk model could predict immunotherapy responses based on three independent cohorts. Patients with low ICD scores showed significantly prolonged survival and were more likely to have clinical benefits from CTLA-4 in Nathanson's cohort ([182]Fig. 7C–E). The low-risk subgroup in the [183]GSE91061 and [184]GSE35640 cohorts exhibited similar therapeutic benefits and clinical responses to ICB therapy ([185]Fig. 7F–I). Furthermore, SubMap analysis confirmed that patients with low ICD scores might exhibit a high likelihood of response to anti-PD1 treatment (Bonferroni corrected P = 0.02, [186]Fig. 7J). Overall, our research strongly suggested that ICD score was a potentially powerful biomarker for predicting the clinical response to immunotherapy. 3.9. Pan-cancer analyses of the ICD score We conducted a series of pan-cancer analyses to elucidate the role of our ICD score across cancers. Generally, the tumors originating from secreting glands (THCA, PCPG, PRAD and ACC) had relatively higher ICD scores, while cancer cells originating from lymph tissue and reproductive organs such as DLBC, OV, UCEC, CESC, and TGCT exhibited lower ones ([187]Fig. 8A and B). Patients were classified into high- and low-risk groups by the median of the population. Subsequent pan-cancer survival analysis demonstrated that samples with high ICD scores in most tumors appeared to be correlated with poor survival indicators. Interestingly, an opposite trend was found in PRAD, UCS, CESC and UCEC. We also estimated some elements associated with tumor immunogenicity including TMB, MSI, SNV Neoantigens, homologous recombination repair deficiency (HRD), chromosomal aneuploidy and intratumor heterogeneity (ITH). Here, we revealed that our ICD score was related to the above parameters in many cancers, such as BRCA, HNSC and LUSC ([188]Figs. S12A–F). Fig. 8. [189]Fig. 8 [190]Open in a new tab Pan-cancer analyses of the ICD score. (A) Box plots (upper panel) displaying the ICD score across cancers ordered using their median values, and bubble chart showing the prognostic value of ICD score in each cancer. (B) Average level of ICD score in individual cancer types. (C) Correlation between ICD score and cancer hallmark signaling pathways across cancers. (D) Correlation between ICD score and the expression of immune regulator genes across cancers. (E) Correlation between ICD score and the level of immune cell infiltration assessed by seven algorithms across cancers. Next, we calculated the enrichment score of hallmark cancer-associated signaling pathways across cancers. The ICD score showed a positive association with immune-inhibited oncogenic pathways ([191]Fig. 8C). Conversely, pronounced negative associations were observed between ICD score and anticancer immune signaling pathways. Similarly, our results depicted a negative correlation between ICD score and the expression of several immune-related genes ([192]Fig. 8D). Seven algorithms (Epic, Quantiseq, Timer, xCell, Mcp-counter, CIBERSORT and CIBERSORT-abs) were used to explore the association between ICD score and immune cell infiltration across pan-cancer. As expected, adverse correlations were observed in the overall TCGA pan-cancer cohort ([193]Fig. 8E). To conclude, our findings revealed that ICD score played a critical role in the immunogenicity and anticancer immunity across multiple cancers. 4. Discussion Despite the significant success of immunotherapy in clinical treatment, the therapeutic effects in HCC remain unsatisfactory and varies greatly among individuals [[194]40,[195]41]. The insufficient tumor immunogenicity and immunosuppressive TME may contribute to the poor clinical outcomes of ICIs [[196]42]. Consequently, there is an urgent need for reliable new predictive biomarkers that can stratify patients and avoid both overtreatment and undertreatment, thereby improving the clinical outcomes of HCC. As a distinctive form of cell demise, ICD can drive the antitumor immunity by releasing TAAs and danger‐associated molecular patterns (DAMPs), ultimately leading to the activation of CTL-driven adaptive immunity and the formation of long-term immune memory [[197]43]. It is reported that treatment-driven ICD reinforces remarkable antitumor effects of conventional anticancer chemotherapies and radiotherapy [[198]44,[199]45]. The combined application of ICD inducers and ICIs shows great potential to reverse immunosuppression, elicit anticancer immune responses and improve the clinical management of HCC patients [[200]46]. Recently, nanoparticles have emerged as a promising vehicle which can target the tumor site and minimize the damage to healthy cells [[201][47], [202][48], [203][49], [204][50]]. Emerging evidence has demonstrated that using nanocarriers for targeted delivery of ICD-inducing drugs to cancer cells can improve the therapeutic outcomes and exhibit synergistic activity with various immunotherapies [[205][51], [206][52], [207][53]]. These studies highlight the important roles of ICD in the tumor immune evasion and alteration of the immune microenvironment. However, the relationship of ICD with immune microenvironment, prognosis, and drug benefits in HCC is not yet fully understood. In this study, we analyzed the multi-omics and clinicopathological data of various cancer types obtained from the TCGA database, revealing that ICD-associated molecules were closely related to clinical outcome and the regulation of genome variation. We identified three distinct ICD subsets characterized by different immune cell infiltration levels, biological behaviors, and prognoses. Cluster A was characterized by the activation of DNA damage repair pathways. Cluster B was correlated with the activation of immune-related pathways and abundant immune infiltration. Cluster C was particularly enriched with unfavorable molecular features such as MYC and TGF-β pathways, with the worst prognosis compared to Cluster A and Cluster B. MYC, being involved in the oncogenic process, is essential for cell proliferation, differentiation, apoptosis and metabolism. The dysregulation of MYC has been reported in up to 70 % of human cancers and is related to adverse survival [[208]54]. Additionally, the multifunctional cytokine TGF-β is widely recognized as an immune-suppressive factor in the TME [[209]55], potentially explaining the poor prognosis observed in this cluster. Considering the pivotal impact of ICD on clinical outcomes and immune regulation in HCC, our study applied a novel computational framework to construct a reliable prediction model. First, the RF algorithm was used to screen out genes with a greater impact on prognosis. Second, to build a more simplified and effective model, we utilized decisive GMM-based clustering, a very feasible method with good clustering performance, to select gene set with the highest AUC value. Afterward, LASSO regularization was applied to create an ICD-related risk score according to the 10 features selected by GMM. Our ICD score was demonstrated to be an independent prognostic factor and patients with low ICD scores exhibited significantly improved survival rates than those with high ICD scores. Its prognostic predictive ability was fully validated in 5 independent HCC cohorts, highlighting the tremendous potential for the clinical application of ICD score. The GSVA results revealed significantly distinct biological functions between two ICD subtypes. We found that a high-risk score was linked to the activation of stromal and carcinogenic pathways, while a low-risk score was correlated with immune activation and metabolism-related pathways. The microenvironment of HCC is a highly complex and dynamic system composed of immune cells, stromal cells, extracellular chemicals, cytokines, and chemokines [[210]56]. Previous studies have confirmed that TME has both positive and negative effects on tumorigenesis, and the magnitude of immune suppression in the TME is closely associated with poor prognosis in HCC patients [[211][57], [212][58], [213][59]]. In our study, patients in low-risk subgroup were found to exhibit significantly higher infiltration levels of anti-tumor immune cells. In line with this, ICD score was negatively linked to the activation of several cancer immunity cycles. Interestingly, the activation of neutrophil recruiting and Treg cell recruiting was observed in high-risk subgroup. Tumor-related neutrophils are crucial to the TME and may trigger tumorous inflammatory response and promote invasion and metastasis [[214]60]. As a member of immunosuppressive cells, Tregs can secrete immunosuppressive cytokines or expressing co-suppressive molecules to inhibit the immune activation [[215]61]. Hence, the recruitment of neutrophil and Treg cell can promote tumor immunosuppression and ultimately lead to a poor immune response. These findings indicated that the low-risk group may represent an “immune-hot” phenotype with abundant immune cell infiltration and powerful anti-tumor immunity. Conversely, patients with high ICD scores may exhibit an “immune-cold” phenotype characterized by low anticancer immunity and bad prognosis. Increasing evidence shows that ICD can convert the TME from immunologically “cold” into “hot”, making them sensitive to immunosuppressive therapy. Hence, we speculate that combined application of ICD inducers and ICIs may have great potential to become an emerging therapeutic strategy and offer promising treatment options for patients in high-risk group. Given the correlation between ICD score and the infiltration characteristics of TME cells, we further investigated the role of ICD score in forecasting the effects of immunotherapy. The TIDE score and submap algorithm revealed that patients with low ICD scores were prone to benefit more from anti-PD-1 treatment. We also validated the predictive value of ICD score in three more independent immunotherapeutic datasets. A significantly difference in ICD scores was observed between non-responders and responders, with the low-risk group showing increased sensitivity to immunotherapy and better survival. RTK inhibitors have been considered as the first-line drug therapies for the treatment of advanced-stage HCC [[216]3]. It was found that ICD score was markedly higher in non-responders compared with responders undergoing sorafenib therapy, indicating that ICD score can effectively discriminate potential patents who may have a better response to RTK inhibitors. Considering the lack of effective treatments and poor prognosis in high-risk group, we identified several potential therapeutic targets and chemotherapy candidates for HCC patients with high ICD scores. Additionally, pan-cancer analysis indicated that ICD score was positively correlated with the activation of immune-inhibited oncogenic pathways and negatively related to the infiltration level of immune cell. These findings further indicated that our ICD score may may be essential for anticancer immunity and act as a valuable generalizable biomarker for evaluating immunotherapy benefits and optimizing the precision treatment, not only in HCC but also across various cancer types. Several limitations in our research need to be addressed. First, although our ICD score was able to predict the prognosis and response to immunotherapy in HCC patients, all samples included in this study were retrospective. Therefore, external validation of the ICD score in a prospective multicenter cohort is needed to determine whether it can replace existing biomarkers. Second, the in-depth molecular mechanisms of how ICD score can affect the TME and immunotherapy response in HCC need further investigation through in vivo and ex vivo research. Third, although we used multiple algorithms to identify promising drug targets and compounds for HCC patients, further experimental and clinical studies are essential to enhance their clinical application. 5. Conclusion In summary, we developed an innovative and reliable model to predict the clinical outcomes, immune infiltration status, and the treatment benefits of chemotherapy, RTK inhibitors, and immunotherapy for HCC patients. Our research offers a roadmap for the stratification of patients, which can optimize the precision therapy and further guide clinical management for individual HCC patients. We believe that our study may provide valuable insights into the role of ICD in HCC and pave the way for new targeted therapies. CRediT authorship contribution statement Dongjing Zhang: Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Formal analysis, Data curation, Conceptualization. Bingyun Lu: Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Data curation. Qianqian Ma: Writing – review & editing, Software, Methodology, Investigation, Data curation. Wen Xu: Software, Methodology, Investigation, Data curation. Qi Zhang: Validation, Investigation, Formal analysis, Data curation. Zhiqi Xiao: Validation, Investigation, Data curation. Yuanheng Li: Validation, Software, Methodology. Ren Chen: Writing – review & editing, Validation, Project administration, Funding acquisition, Conceptualization. An-jiang Wang: Writing – review & editing, Validation, Project administration, Funding acquisition, Conceptualization. Ethics declarations Informed consent was not required for this study. Data availability statement All data underlying this article can be downloaded from The Cancer Genome Atlas ([217]https://portal.gdc.cancer.gov/), Gene Expression Omnibus (https:// [218]www.ncbi.nlm.nih.gov/geo/), and International Cancer Genome Consortium (https:// dcc.icgc.org/). Additional information can be requested from the corresponding authors. Funding This study was supported by the National Natural Science Foundation of China (No. 82160115), Shenzhen Medical Research Funding (No. A2301050), Jiangxi Provincial Science and Technology Department key research and development project (20202BBGL73093), Shenzhen Clinical Research Center for Digestive Disease, Shenzhen Hospital, Southern Medical University, Shenzhen, China. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements