
==== Front
NPJ Precis Oncol
NPJ Precis Oncol
NPJ Precision Oncology
2397-768X
Nature Publishing Group UK London

39245753
693
10.1038/s41698-024-00693-9
Article
Immune-related cell death index and its application for hepatocellular carcinoma
Sun Zhao 1
Liu Hao 1
Zhao Qian 2
Li Jie-Han 1
Peng San-Fei 1
Zhang Zhen 1
Yang Jing-Hua jhy@zzu.edu.cn

2
Fu Yang fuyang@zzu.edu.cn

1
1 https://ror.org/056swr059 grid.412633.1 Department of Gastrointestinal Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, Henan China
2 https://ror.org/056swr059 grid.412633.1 Clinical Systems Biology Key Laboratory, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, Henan China
8 9 2024
8 9 2024
2024
8 19424 11 2023
28 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Regulated cell death (RCD) plays a crucial role in the immune microenvironment, development, and progression of hepatocellular carcinoma (HCC). However, reliable immune-related cell death signatures have not been explored. In this study, we collected 12 RCD modes (e.g., apoptosis, ferroptosis, and cuproptosis), including 1078 regulators, to identify immune-related cell death genes based on HCC immune subgroups. Using a developed competitive machine learning framework, nine genes were screened to construct the immune-related cell death index (IRCDI), which is available for online application. Multi-omics data, along with clinical features, were analyzed to explore the HCC malignant heterogeneity. To validate the efficacy of this model, more than 18 independent cohorts, including survival and diverse treatment cohorts and datasets, were utilized. These findings were further validated using in-house samples and molecular biological experiments. Overall, the IRCDI may have a wide application in individual therapeutic decision-making and improving outcomes for HCC patients.

Subject terms

Cancer models
Cancer microenvironment
https://doi.org/10.13039/501100001809 National Natural Science Foundation of China (National Science Foundation of China) 81871995 Fu Yang issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Globally, hepatocellular carcinoma (HCC) ranks third in cancer-related deaths due to late-stage diagnosis and limited treatment options1,2. Targeted drugs and other treatments, including immune checkpoint inhibitors (ICIs) and trans-arterial chemoembolization (TACE), have shown great potential in HCC treatment3,4. However, the predictive efficacy of these therapies has consistently been unsatisfactory4. In particular, predicting patient response to immunotherapy is challenging due to the complexity of the immune response, which includes the status of immune checkpoints, intercellular communication, immune cell infiltration, intrinsic mechanisms of tumor cells, and the tumor microenvironment (TME)3–7. Therefore, there is an urgent need to develop reliable biomarkers and models that can accurately distinguish HCC subtypes and identify patients who are most likely to benefit from drug-based therapies. Exploration of tumor features by generating multi-biological processes can help overcome this challenge.

Currently, a total of 12 regulated cell death (RCD)/programmed cell death (PCD) modes have been identified in cell cycle regulation, and accumulating evidence suggests that cell death plays a key role in cancer immunity8. Autophagy is a tightly regulated cell death pathway characterized by stress-induced catabolism. During autophagy, adenosine triphosphate (ATP) is secreted into the extracellular microenvironment, mediating robust anti-cancer immune responses9. Pyroptosis, executed by the gasdermin (GSDM) family, is a lytic and inflammatory type of programmed cell death. Dying cancer cells can release a large number of damage-associated molecular patterns (DAMPs) during pyroptosis, robustly initiating anti-cancer immune responses10. Ferroptosis, a type of immune cell death (ICD), also involves the release of multiple DAMPs in the early stage, inducing adaptive immune responses11. Apart from these well-known RCD modes, autophagy and apoptosis induced by tumor necrosis factor (TNF) evoke anti-cancer immune responses, revealing novel strategies for immunotherapy12. Current studies have shown the positive feedback loop between cell death and immune response13,14. Consequently, some approved drugs (e.g., cell death inducers: PEM/CDDP/MEKi) combined with immune checkpoint inhibitors (ICIs) have been clinically applied, demonstrating remarkable effects such as direct tumor shrinkage and promoting immune cell infiltration in TME15. Although tight associations have been identified, other cell death modes, including cuproptosis, NETosis, alkaliptosis, parthanatos, lysosome-dependent cell death (LCD), and oxeiptosis, are still poorly characterized16.

Here, to decipher the tumor heterogeneity and improve individual therapeutic assessment in HCC, we investigated 12 regulated cell death modes, including 1078 regulators, and concluded a gene collection termed immune-related cell death genes (IRCDGs) based on our defined immune-cold- and hot- HCC subgroups. Subsequently, a competitive machine learning framework was employed to construct an immune-related cell death index (IRCDI) consisting of nine genes selected from the IRCDGs. These key genes were validated in our cohort and subjected to deeper analysis using proteomics. Among them, ACSL6, a member of the ACSL family protein associated with ferroptosis17, was confirmed to act as a protective protein in HepG2 and Huh7 cell lines through laboratory experiments and mass spectrometry technology. The IRCDI was further validated using more than 18 independent databases to predict outcomes, immune responses, and drug sensitivity in HCC patients.

Results

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).Fig. 1 Identification and validation of two distinct immune subgroups in HCC.

a Clustering heatmap showing consensus scores of HCC samples with two subgroups. b Patients were classified by consensus and M3C cluster method. c CDF curves with all numbers of subgroups (k = 2–9). d The proportion of ambiguous clustering (PAC) score indicating the optimal number k is 2. e Heatmap showing the abundance of 28 immune cells calculated by ssGSEA method in two distinct immune subtypes. f–i Heatmaps and PCA plots showing the IRCDGs expression pattern and distributions of immune subgroups in independent cohorts. j The cold- and hot-immune subtypes in four published studies in the TCGA-LIHC cohort.

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 studies18–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 pathway18,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).Fig. 2 Genomics and functional features of immune-related cell death genes (IRCDGs).

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 (https://theshy.shinyapps.io/IRCDI/) 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).Fig. 3 Construction of immune-related cell death index (IRDCI) by developing the competitive machine learning framework.

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).Fig. 4 Clinical features associations and prognosis prediction of IRCDI.

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 < 0.05; **P < 0.01; ***P < 0.001.

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).Fig. 5 Distinct driver genes and biological pathways in IRCDI subgroups.

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 < 0.0001, Fig. 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 < 0.0001). The green module showed a significantly negative correlation to IRCDI (p-value < 0.0001). The module memberships of the two gene modules were displayed in Supplementary Fig. 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.Fig. 6 Immune cell infiltration and PD-1 associations with IRCDI.

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 < 0.05; **P < 0.01; ***P < 0.001. 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 < 0.0001, Supplementary Fig. 12)23. 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 < 0.0001, Supplementary Fig. 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 cells24. 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.Fig. 7 Treatment response performance of IRCDI.

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.Fig. 8 FDA-approved common drug selection according to IRCDI.

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.

Fig. 9 Experimental exploration of the nine-panel genes.

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 < 0.05; **P < 0.01; ***P < 0.001. Centre line: the median; Bounds of the box: the 25th and 75th percentile.

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.

Discussion

HCC poses a serious threat to global public health1. For advanced HCC, multiple treatment modalities are now available, including percutaneous ablation, radiation, liver transplantation, and target drug therapy. Among them, ICIs have developed rapidly as a drug therapy in recent years25. However, due to the high degree of tumor heterogeneity, there is great variation in the response of HCC to treatment26. Activation of different modes of cell death, which underlie tumor shrinkage, is key to immunotherapy. Pyroptosis, apoptosis, and necroptosis are previously regarded as forms of inflammatory programmed cell death associated with innate immunity, and specific products of dead cells serve as important drivers of immune response, inflammation, and autoreactivity27. Later studies summarize this phenomenon, termed immune cell death (ICD), as a “form of regulated cell death which is sufficient to activate an adaptive immune response in immunocompetent syngeneic hosts” and emphasized the importance of ICD inducers in tumor therapy28,29. However, owing to the streetlight effect, current studies focus on the typically regulated cell death modes and thus lack the landscape of the whole picture between cell death and the immune microenvironment. To solve this problem, some researchers have started to apply the combined analysis of multi-defined biological processes, such as cuproptosis and ferroptosis, which share the consensus pathways, including glutathione metabolism and oxidative stress injury functioning in the regulation of the metal ions30. Here, to cover the less-charactered cell death modes, we systemically analyzed 12 modes of regulated cell death involving 1078 genes in HCC immune subtypes to identify the signature that predicts drug benefit for improved prognosis.

For an unbiased analysis of 12 cell death modes, we constructed an immune-related cell death index (IRCDI) by a competitive machine learning framework based on the IRCDGs. This strategy showed high interpretability and robustness and has been proven effective in previous studies31. A total of nine model genes, including six risk genes (PPT1, HSPB8, HMOX1, NQO1, SERPINE1, and FGB) and three protective genes (ACSL6, FYN, and LGALS9), were identified as core IRCDG and validated by qPT-PCR by in-house cohort. Among them, several genes are reported to function in regulating cell death and immune response. Palmitoyl-protein thioesterase 1 (PPT1) is considered an unfavorable gene in HCC and is up-regulated in sorafenib-resistant cell lines. The effect related to sorafenib sensitivity was then confirmed by DC661 (PPT1 inhibitor) which inhibits autophagy and induces mitochondrial pathway apoptosis32. Sharma et al. also reported that PPT1 inhibition can enhance the efficacy of anti-PD-1 drugs33, which is consistent with our results. Quinone oxidoreductase 1 (NQO1) plays an important role in proper innate sensing of the tumor microenvironment, limiting the efficacy of immunotherapies targeting T cells34. The combination of NQO1 inhibitors and immunotherapy can have a huge effect on overcoming adaptive drug resistance34,35. In breast cancer, depletion of heme oxygenase 1 (HMOX1), an immunosuppressive gene, induces T cell-mediated cancer cytotoxicity by eliminating the ability of cancer cells to activate Tregs and overcome in situ immune attack36. Patients with high SERPINE1 expression tend to exhibit poorer prognosis and lower immune cell infiltration in renal clear cell carcinoma37. Of note, ACSL6 was identified as a protective protein in HCC. Knocking down the ACSL6 changed pathways, including the well-known markers such as L1CAM and RFTN1, which functioned in tumor proliferation and metastasis38. Overall, all model genes from cell death regulators demonstrated a strong association between cell death and immune response, suggesting potential drug targets in HCC.

American Joint Committee on Cancer (AJCC) staging system is used to evaluate outcomes in cancer patients and provide guidelines for intervention39. However, this system cannot meet additional clinical requirements, such as drug response prediction. Machine learning has been applied to predict therapy responses in colorectal and breast cancer31,40. Our developed IRCDI showed better performance than the AJCC system in prognostic prediction. Furthermore, the IRCDI was negatively correlated with immune checkpoint PD-1 in multiple independent cohorts, which was validated in our HCC samples. Thus, we investigated its ability to predict immunotherapy. The TIDE and four independent immunotherapy cohorts’ results demonstrated that the model achieved satisfactory accuracy in drug response prediction. Furthermore, this model performed well even in predicting the TACE treatment response. All the results showed the IRCDI had a strong generalization ability, owing to the wide involvement of immunity and cell death processes, which were directly relevant to any treatment modality16,41,42. Moreover, the model was used to screen possible antitumor drugs in HCC. After evaluating 198 common FDA-approved drugs, sepantronium bromide (YM155) and AZD7762 were identified as potential drug targets by drug response prediction and experimental recorded database (GDSC), which showed antitumor effects in other cancers43,44. The high correlation between IRCDI and target gene expression (BIRC5 and CHK1) supports the drug screening results.

Our study has several limitations. Firstly, although our results were based on more than 18 multi-independent datasets, all included samples were from retrospective survival datasets, and thus this model should be validated by more prospective cohorts. Secondly, we evaluated the relationship between cell death and immune cell infiltration based on regulators, which may result in bias for less-characterized cell death types. Lastly, in vivo research is needed for the exploration of the model genes’ biological functions.

In conclusion, we revealed associations between diverse modes of cell death and immune subtypes and constructed a robust model (IRCDI) to predict HCC patient prognosis and sensitivity to diverse treatments. We released this model online (https://theshy.shinyapps.io/IRCDI/), thereby providing a useful tool for clinicians.

Methods

Datasets collection and processing

A total of 12 cell death modes were analyzed involving 1078 regulators obtained from the MSigDB database and published articles after passing manual review45,46. All genes could be retrieved in Supplementary Data 1. HCC and other tumor transcriptome profiles, genomic data, and clinical data were retrieved from the Cancer Genome Atlas (TCGA) database, the International Cancer Genome Consortium (ICGC) database, cBioPortal, and the Gene Expression Omnibus (GEO) databases using the TCGAbiolink, cBioProtal, and GEOquery R packages. Overall, 19 non-redundant datasets were analyzed, consisting of 8 survival datasets (TCGA-LIHC, ICGC-LIR1, GSE54236, GSE14520, TCGA-PAAD, TCGA-STAD, TCGA-ESCA, and TCGA-COAD), one large transcriptome dataset (GSE25097), four immunotherapy datasets (GSE100797, IMvigor 210 cohort, GSE35640, and GSE91061), one TACE treatment dataset (GSE104580), one proteomics dataset (In-house HCC cohort), one single-cell dataset GSE149614, and 3 HCC cell lines datasets. The HCC cell lines transcriptomics and the IC50 of the common drugs were downloaded from the Cancer Cell Line Encyclopedia (CCLE) database and the Genomics of Drug Sensitivity in Cancer (GDSC) database. The HCC line lines proteomics data were obtained from our experiments. Among these cohorts, TCGA-LIHC was treated as the training dataset, and other datasets were regarded as the validating datasets. All details, including number of patients/cells and platform, could be retrieved in Supplementary Data 2.

In this study, RNA expression profiles of included datasets from next-generation sequencing (Illumina platform based on NGS technology, Supplementary Data 2) were transformed into transcripts per million (TPM) and log-2 transformed. RNA expression profiles by array were processed by the robust multiarray averaging (RMA) method using the “Affy” package. The gene annotation of probes from different platforms used the “idmap3” package. Before constructing the model, all gene expression data were scaled to Z-score value before model construction via the “scale” function in the R “base” package, which maintained most gene expression levels within a ±2 range. Data from different cohorts were analyzed separately to construct and validate the model.

Immune cell filtration evaluation

To measure the relative abundance of immune cells, single-sample gene set enrichment analysis (ssGSEA) was applied using the R package GSVA via the “gsva” function47. To validate the robustness of ssGSEA immune cell infiltration results, seven other algorithms were used, including ESTIMATE, quanTIseq, xCell, MCP-Counter, CIBERSORT, CIBERSORT-ABS, and TIMER48–53. Detail parameters were described in the main parameters section.

Consensus clustering, validation, and comparison subgroups with published studies

We applied the consensus clustering algorithm in the R package consensusClusterPlus via “ConsensusClusterPlus” function to identify distinct immune subgroups in HCC54. The optimum subgroup number was simultaneously defined by using cumulative distribution function curves, the method of proportion of ambiguous clustering (PAC), together with the “Best.nc” of NbClust chosen by diverse criteria55. To validate the robustness of the immune subtypes, M3C (Monte Carlo Reference-based Consensus Clustering) clustering by M3C package using the default parameters56. Besides, the relative cluster stability index (RCSI) was used to confirm selecting of the optimal number of subgroups. Detail parameters were described in the main parameters section.

Four published studies were used to compare our immune subtype. According to the uploaded gene set from these studies, we reproduced the HCC clusters defined by these studies in the TCGA-LIHC cohort as far as possible. Then, we calculated the proportions of immune subgroups in each cluster in each study. The feature genes could be obtained from Supplementary Data 4.

To validate the structural robustness of the immune subgroups, hierarchical clustering was first performed to divide the dataset into two subgroups based on the IRCDGs, which were regarded as the immune signature genes in this study (Supplementary Data 3). Then SubMap analysis (an unsupervised subclass mapping method) was performed to calculate the statistical significance in the TCGA-LIHC cohort and three other cohorts.

Expression and variation levels of immune-related cell death genes

Cell death genes differentially expressed in distinct immune subtypes were regarded as immune-related cell death genes (IRCDGs) (significance level: adjusted p < 0.05 and |log2FC | > 0.5) using the “limma” R package. Chromosomal locations of IRCDGs were labeled using the “circlize” package. Somatic mutations of IRCDGs in HCC were visualized using the “Maftools” package57.

Pathway activation evaluation and protein-protein interaction analysis

Different genes were annotated using the “clusterProfiler” package in R via the “enrichGO” and “enrichKEGG” functions to explore the enrichment pathways58. The protein-protein interaction analysis (PPI) was performed by Metascape using default parameters59. Gene set enrichment analysis (GSEA) via the “GSEA” function using default parameters was applied to determine differential functions in low- and high-IRCDI groups. The results were visualized via “gseaplot2” function in the “enrichplot” package. The tumor hallmark activation evaluation was calculated using the ssGSEA methods described. The gene collections were obtained from MsigDB (https://www.gsea-msigdb.org/gsea/msigdb). Then, we performed Spearman correlation analysis between quantified pathways and IRCDI to explore the most relevant pathways.

Construction of the immune-related cell death index

We constructed a prognostic model of IRCDGs by a competitive machine learning framework, with the final score of each patient regarded as the IRCDI. Step 1: Univariate Cox regression: Hazard ratios (HRs) with the 95% confidence intervals (CIs) of each IRCDG were first determined by univariate Cox regression analysis, of which, genes with p-value less than 0.05 were regarded as prognostic genes. Step 2: Least absolute shrinkage and selection operator (LASSO) Cox regression: To avoid overfitting, LASSO Cox regression was used to reduce the number of candidate genes using the “glmnet” package. Those genes that remained in the model had minimum lambda values. Step 3: Competitive machine learning framework. This framework consisted of nine models including 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), which returned the predicted score and then calculated the Harrell’s concordance index (C-index) of each cohort. The detailed parameters of these machine learning models could be obtained from the data availability section. For StepCox method: containing survival genes with minimum Akaike information criterion (AIC) was regarded as the final model and exported into the IRCDI using the following formula (1):1 IRCDI=∑i=1nCoefi*Expi

Where N denotes the number of model genes, Coefi denotes each model gene risk coefficient, and Expi denotes each model gene mRNA level.

Weighted correlation network analysis

We employed weighted correlation network analysis (WGCNA) to explore the co-expression genes of IRCDI. An optimal soft threshold β was determined to satisfy the requirements for a scale-free network. Subsequently, the weighted adjacency matrix was transformed into a topological overlap matrix (TOM), and the corresponding dissimilarity (1- TOM) was computed. The dynamic tree-cutting approach was utilized for module identification. To identify gene modules significantly associated with IRCDI, we selected the module with the highest correlation for further biological pathway enrichment investigation.

Evaluation of IRCDI and establishment of the online nomogram

Harrell’s concordance index (C-index) and 95% CI of each dataset were calculated by the “concordance.index” function using the default parameters in “pec” R package. Time-dependent areas under the curve (AUC) were applied to internal and external datasets using the “timeROC” function in the “timeROC” package. The 95% CIs of the AUCs were calculated by the “confint” function (level = 0.95) in the same package. Nomogram and calibration curves were calculated and displayed using the “rms” and “ggplot2” packages. The dynamic nomogram was constructed using the “shinyPredict” R package and then deployed to the server as a web application. Detail parameters were described in the main parameters section.

Therapeutic response prediction

Predicted half maximal inhibitory concentrations (IC50) of 198 common antitumor drugs in HCC were calculated using the oncoPredict R package60. The Tumor Immune Dysfunction and Exclusion (TIDE) score was used to predict responses of ICIs (https://tide.dfci.harvard.edu/).

Cell culture, treatment, CCK8, and wound healing assay

Human hepatocellular carcinoma cell lines HepG2 and Huh7 were cultured in the recommended growth medium supplemented with 10% fetal bovine serum (FBS). The cells were maintained in a 5% CO2 atmosphere at 37 °C. Transfection of cells with siRNA was performed using Lipofectamine 2000 (Life Technologies, USA) according to the manufacturer’s instructions. Nonspecific siRNAs were used as negative controls. The siRNAs were designed and synthesized by Sangon Biotech (China).

HepG2 and Huh7 cells with respective treatments were seeded at a density of 2*103 cells/well in 96-well cell culture plates. Cell viability was assessed using the CCK-8 assay at 1-, 2-, 3-, and 4 days post-treatment, following the manufacturer’s instructions. The absorbance of the cells at 450 nm (OD450) was measured using a 96-well plate reader (DG5032, Hua Dong, Nanjing, China).

For the wound healing assay, transfected HepG2 and HuH7 cells were seeded in 6-well plates at a density of 1*105 cells/well. Once the cells reached over 90% confluence, a wound was created by scratching with a 200 μL plastic pipette tip. The width of the wound was recorded and photographed at 0 and 48 hours.

Quantification real-time reverse transcriptase (qRT-PCR)

The HCC tissues were collected in the First Affiliated Hospital of Zhengzhou University (China), which was approved by the ethics committee of this hospital (2021-KY-0852-004). Total RNA was extracted from tissues and cell lines using Tri-Reagent (TIANGEN BIOTECH, Beijing, China) based on the manufacturer’s protocol. The quantity and the RNA purity were evaluated using spectrophotometry, and 3 μg of RNA from each sample was used for the RT reaction. Transcripts of β-actin were used as internal controls for normalization purposes.

LC–MS/MS analysis

All tissues and cell lines were lysed in RIPA buffer (Merck) supplemented with protease and phosphatase inhibitors. The lysates were centrifuged, and the supernatants were collected. The concentration of protein lysates was measured by the BCA Protein Assay Kits (Beyotime). Equal amounts of the lysates were precipitated with pre-cooled acetone at −80 °C. The protein precipitates were resuspended in 25 mM ammonium bicarbonate solution. Proteins were digested with trypsin. Afterward, peptides were desalted by Sep-Pak C18 Vac cartridges (Waters) and freeze-dried in a Speed Vac device (Thermo Fisher Scientific). Finally, the peptides were re-dissolved in 0.1% formic acid water and quantified by Pierce Quantitative Colorimetric Peptide Assay Kit.

The LC–MS/MS system consisted of a nanoflow high-performance liquid chromatograph instrument (EASY-nLC1200 nanofowLC, Thermo Fisher Scientific) coupled to a thermo Q Exactive HF-X-Orbitrip mass spectrometer (Thermo Fisher Scientific). For data acquisition, peptides were separated on the analytical column (75 μm × 150 mm, C18-AQ, 1.9 μm), and the total elution time was 120 min. Then, the raw files were searched against the human Uniprot database (https://www.uniprot.org/) with Protein Discovery (version 2.4). The digestion mode was set to specific, and trypsin/P was chosen. Oxidation of methionine and acetylation of any N-terminal were set as variable modifications. The false discovery rate (FDR) value was set to 1%.

Main parameters

For ssGSEA analysis, the “gsva” function parameters were set as mx. diff = T, method = “ssgsea”, and kcdf = “Gaussian”. For ConsensusClusterPlus analysis, the parameters were set as reps = 1000, pItem = 0.8, pFeature = 1, distance = “euclidean”, clusterAlg = “km”, and maxK = 9. For timeROC analysis, the parameters were set as weighting = ‘marginal’, ROC = T, and iid = T.

Statistical analysis

All statistical analyses were conducted using R software (v4.1.0). The continuous variables were reported as the Mean ± SD, which had been widely shown in the box plots. Limma, Wilcoxon, and T-tests were applied to analyze differences between the two groups according to the different data types. Kruskal–Wallis test was applied to analyze differences between multiple groups. Meanwhile, a paired T-test was used to compare the data from cancer and paired adjacent tissue. Cut-off points of the continuous variables (IRCDI and genes) were determined by the “surv_cutpoint” function in the “survminer” package and the log-rank test was used to determine survival differences using the “survival” package. Except for survival analysis, the HCC were divided into low- and high-IRCDI groups by the median of IRCDI.

In this study, all correlation analyses use Spearman correlation via the “cor.test” function in the R base package. Two-sided p-values < 0.05 were considered statistically significant. The adjusted p-values from the differential expression analysis and pathway enrichment analysis were used by the Benjamin Hochberg (BH) method.

Supplementary information

Supplementary information

Supplementary Data 1

Supplementary Data 2

Supplementary Data 3

Supplementary Data 4

Supplementary information

The online version contains supplementary material available at 10.1038/s41698-024-00693-9.

Acknowledgements

This study was supported by the Scientific Research and Innovation Team of the First Affiliated Hospital of Zhengzhou University (QNCXTD2023022), the Excellent Youth Fund of Henan Natural Science Foundation (212300410075), the Project of Henan Health Department (YXKC2020029), and the Chinese National Natural Science Foundation (Nos. 81871995).

Author contributions

Y.F. and J.Y. contributed to the conception and design of the study. Z.S. organized the database. Z.S. and H.L. performed the statistical analysis. H.L., Q.Z., and J.L. applied lab experiments. Z.S., H.L., and Q.Z. wrote the first draft of the manuscript. S.P., Z.Z., Y.F., and J.Y. wrote sections of the manuscript. All authors contributed to the manuscript revision, read, and approved the submitted version.

Data availability

All datasets in this study can be obtained from the Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and International Cancer Genome Consortium (ICGC) database. The essential R scripts can be obtained from GitHub (https://github.com/FuLab-ZhaoSun/IRCDI/). All other data, including proteomics data of in-house HCC, can be acquired from the corresponding authors upon reasonable request. To facilitate researchers, this model has been made available online (https://theshy.shinyapps.io/IRCDI/).

Competing interests

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Zhao Sun, Hao Liu, Qian Zhao.
==== Refs
References

1. Vogel A Hepatocellular carcinoma Lancet 2022 400 1345 1362 10.1016/S0140-6736(22)01200-4 36084663
Vogel, A. et al. Hepatocellular carcinoma. Lancet 400, 1345–1362 (2022).36084663 10.1016/S0140-6736(22)01200-4
2. Konyn P Ahmed A Kim D Current epidemiology in hepatocellular carcinoma Expert Rev. Gastroenterol. Hepatol. 2021 15 1295 1307 10.1080/17474124.2021.1991792 34624198
Konyn, P., Ahmed, A. & Kim, D. Current epidemiology in hepatocellular carcinoma. Expert Rev. Gastroenterol. Hepatol. 15, 1295–1307 (2021).34624198 10.1080/17474124.2021.1991792
3. Huang A Targeted therapy for hepatocellular carcinoma Signal Transduct. Target Ther. 2020 5 146 10.1038/s41392-020-00264-x 32782275
Huang, A. et al. Targeted therapy for hepatocellular carcinoma. Signal Transduct. Target Ther. 5, 146 (2020).32782275 10.1038/s41392-020-00264-x
4. Yang C Evolving therapeutic landscape of advanced hepatocellular carcinoma Nat. Rev. Gastroenterol. Hepatol. 2023 20 203 222 10.1038/s41575-022-00704-9 36369487
Yang, C. et al. Evolving therapeutic landscape of advanced hepatocellular carcinoma. Nat. Rev. Gastroenterol. Hepatol. 20, 203–222 (2023).36369487 10.1038/s41575-022-00704-9
5. Craig AJ Tumour evolution in hepatocellular carcinoma Nat. Rev. Gastroenterol. Hepatol. 2020 17 139 152 10.1038/s41575-019-0229-4 31792430
Craig, A. J. et al. Tumour evolution in hepatocellular carcinoma. Nat. Rev. Gastroenterol. Hepatol. 17, 139–152 (2020).31792430 10.1038/s41575-019-0229-4
6. Li L Wang H Heterogeneity of liver cancer and personalized therapy Cancer Lett. 2016 379 191 197 10.1016/j.canlet.2015.07.018 26213370
Li, L. & Wang, H. Heterogeneity of liver cancer and personalized therapy. Cancer Lett. 379, 191–197 (2016).26213370 10.1016/j.canlet.2015.07.018
7. Oura K Tumor immune microenvironment and immunosuppressive therapy in hepatocellular carcinoma: a review Int. J. Mol. Sci. 2021 22 5801 10.3390/ijms22115801 34071550
Oura, K. et al. Tumor immune microenvironment and immunosuppressive therapy in hepatocellular carcinoma: a review. Int. J. Mol. Sci. 22, 5801 (2021).34071550 10.3390/ijms22115801
8. Galluzzi L Immunogenic cell death in cancer and infectious disease Nat. Rev. Immunol. 2017 17 97 111 10.1038/nri.2016.107 27748397
Galluzzi, L. et al. Immunogenic cell death in cancer and infectious disease. Nat. Rev. Immunol. 17, 97–111 (2017).27748397 10.1038/nri.2016.107
9. Michaud M Autophagy-dependent anticancer immune responses induced by chemotherapeutic agents in mice Science 2011 334 1573 1577 10.1126/science.1208347 22174255
Michaud, M. et al. Autophagy-dependent anticancer immune responses induced by chemotherapeutic agents in mice. Science 334, 1573–1577 (2011).22174255 10.1126/science.1208347
10. Rosenbaum SR Wilski NA Aplin AE Fueling the fire: inflammatory forms of cell death and implications for cancer immunotherapy Cancer Discov. 2021 11 266 281 10.1158/2159-8290.CD-20-0805 33451983
Rosenbaum, S. R., Wilski, N. A. & Aplin, A. E. Fueling the fire: inflammatory forms of cell death and implications for cancer immunotherapy. Cancer Discov. 11, 266–281 (2021).33451983 10.1158/2159-8290.CD-20-0805
11. Efimova I Vaccination with early ferroptotic cancer cells induces efficient antitumor immunity J. Immunother. Cancer 2020 8 e001369 10.1136/jitc-2020-001369 33188036
Efimova, I. et al. Vaccination with early ferroptotic cancer cells induces efficient antitumor immunity. J. Immunother. Cancer 8, e001369 (2020).33188036 10.1136/jitc-2020-001369
12. Annibaldi A Meier P Checkpoints in TNF-induced cell death: implications in inflammation and cancer Trends Mol. Med. 2018 24 49 65 10.1016/j.molmed.2017.11.002 29217118
Annibaldi, A. & Meier, P. Checkpoints in TNF-induced cell death: implications in inflammation and cancer. Trends Mol. Med. 24, 49–65 (2018).29217118 10.1016/j.molmed.2017.11.002
13. Liao P CD8(+) T cells and fatty acids orchestrate tumor ferroptosis and immunity via ACSL4 Cancer Cell 2022 40 365 378.e6 10.1016/j.ccell.2022.02.003 35216678
Liao, P. et al. CD8(+) T cells and fatty acids orchestrate tumor ferroptosis and immunity via ACSL4. Cancer Cell 40, 365–378.e6 (2022).35216678 10.1016/j.ccell.2022.02.003
14. Zhang Z Gasdermin E suppresses tumour growth by activating anti-tumour immunity Nature 2020 579 415 420 10.1038/s41586-020-2071-9 32188940
Zhang, Z. et al. Gasdermin E suppresses tumour growth by activating anti-tumour immunity. Nature 579, 415–420 (2020).32188940 10.1038/s41586-020-2071-9
15. Limagne E MEK inhibition overcomes chemoimmunotherapy resistance by inducing CXCL10 in cancer cells Cancer Cell 2022 40 136 152.e12 10.1016/j.ccell.2021.12.009 35051357
Limagne, E. et al. MEK inhibition overcomes chemoimmunotherapy resistance by inducing CXCL10 in cancer cells. Cancer Cell 40, 136–152.e12 (2022).35051357 10.1016/j.ccell.2021.12.009
16. Gao W Autophagy, ferroptosis, pyroptosis, and necroptosis in tumor immunotherapy Signal Transduct. Target Ther. 2022 7 196 10.1038/s41392-022-01046-3 35725836
Gao, W. et al. Autophagy, ferroptosis, pyroptosis, and necroptosis in tumor immunotherapy. Signal Transduct. Target Ther. 7, 196 (2022).35725836 10.1038/s41392-022-01046-3
17. Quan J Bode AM Luo X ACSL family: the regulatory mechanisms and therapeutic implications in cancer Eur. J. Pharmacol. 2021 909 174397 10.1016/j.ejphar.2021.174397 34332918
Quan, J., Bode, A. M. & Luo, X. ACSL family: the regulatory mechanisms and therapeutic implications in cancer. Eur. J. Pharmacol. 909, 174397 (2021).34332918 10.1016/j.ejphar.2021.174397
18. Boyault S Transcriptome classification of HCC is related to gene alterations and to new therapeutic targets Hepatology 2007 45 42 52 10.1002/hep.21467 17187432
Boyault, S. et al. Transcriptome classification of HCC is related to gene alterations and to new therapeutic targets. Hepatology 45, 42–52 (2007).17187432 10.1002/hep.21467
19. Chiang DY Focal gains of VEGFA and molecular classification of hepatocellular carcinoma Cancer Res. 2008 68 6779 6788 10.1158/0008-5472.CAN-08-0742 18701503
Chiang, D. Y. et al. Focal gains of VEGFA and molecular classification of hepatocellular carcinoma. Cancer Res. 68, 6779–6788 (2008).18701503 10.1158/0008-5472.CAN-08-0742
20. Hoshida Y Molecular classification and novel targets in hepatocellular carcinoma: recent advancements Semin. Liver Dis. 2010 30 35 51 10.1055/s-0030-1247131 20175032
Hoshida, Y. et al. Molecular classification and novel targets in hepatocellular carcinoma: recent advancements. Semin. Liver Dis. 30, 35–51 (2010).20175032 10.1055/s-0030-1247131
21. Coulouarn C Factor VM Thorgeirsson SS Transforming growth factor-beta gene expression signature in mouse hepatocytes predicts clinical outcome in human cancer Hepatology 2008 47 2059 2067 10.1002/hep.22283 18506891
Coulouarn, C., Factor, V. M. & Thorgeirsson, S. S. Transforming growth factor-beta gene expression signature in mouse hepatocytes predicts clinical outcome in human cancer. Hepatology 47, 2059–2067 (2008).18506891 10.1002/hep.22283
22. Malladi S Metastatic latency and immune evasion through autocrine inhibition of WNT Cell 2016 165 45 60 10.1016/j.cell.2016.02.025 27015306
Malladi, S. et al. Metastatic latency and immune evasion through autocrine inhibition of WNT. Cell 165, 45–60 (2016).27015306 10.1016/j.cell.2016.02.025
23. Hegde S Leader AM Merad M MDSC: markers, development, states, and unaddressed complexity Immunity 2021 54 875 884 10.1016/j.immuni.2021.04.004 33979585
Hegde, S., Leader, A. M. & Merad, M. MDSC: markers, development, states, and unaddressed complexity. Immunity 54, 875–884 (2021).33979585 10.1016/j.immuni.2021.04.004
24. Pinato DJ Trans-arterial chemoembolization as a loco-regional inducer of immunogenic cell death in hepatocellular carcinoma: implications for immunotherapy J. Immunother. Cancer 2021 9 e003311 10.1136/jitc-2021-003311 34593621
Pinato, D. J. et al. Trans-arterial chemoembolization as a loco-regional inducer of immunogenic cell death in hepatocellular carcinoma: implications for immunotherapy. J. Immunother. Cancer 9, e003311 (2021).34593621 10.1136/jitc-2021-003311
25. Llovet JM Locoregional therapies in the era of molecular and immune treatments for hepatocellular carcinoma Nat. Rev. Gastroenterol. Hepatol. 2021 18 293 313 10.1038/s41575-020-00395-0 33510460
Llovet, J. M. et al. Locoregional therapies in the era of molecular and immune treatments for hepatocellular carcinoma. Nat. Rev. Gastroenterol. Hepatol. 18, 293–313 (2021).33510460 10.1038/s41575-020-00395-0
26. Ng HHM Immunohistochemical scoring of CD38 in the tumor microenvironment predicts responsiveness to anti-PD-1/PD-L1 immunotherapy in hepatocellular carcinoma J. Immunother. Cancer 2020 8 e000987 10.1136/jitc-2020-000987 32847986
Ng, H. H. M. et al. Immunohistochemical scoring of CD38 in the tumor microenvironment predicts responsiveness to anti-PD-1/PD-L1 immunotherapy in hepatocellular carcinoma. J. Immunother. Cancer 8, e000987 (2020).32847986 10.1136/jitc-2020-000987
27. Christgen S Tweedell RE Kanneganti TD Programming inflammatory cell death for therapy Pharmacol. Ther. 2022 232 108010 10.1016/j.pharmthera.2021.108010 34619283
Christgen, S., Tweedell, R. E. & Kanneganti, T. D. Programming inflammatory cell death for therapy. Pharmacol. Ther. 232, 108010 (2022).34619283 10.1016/j.pharmthera.2021.108010
28. Galluzzi L Consensus guidelines for the definition, detection and interpretation of immunogenic cell death J. Immunother. Cancer 2020 8 e000337 10.1136/jitc-2019-000337 32209603
Galluzzi, L. et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J. Immunother. Cancer 8, e000337 (2020).32209603 10.1136/jitc-2019-000337
29. Galluzzi L Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018 Cell Death Differ. 2018 25 486 541 10.1038/s41418-017-0012-4 29362479
Galluzzi, L. et al. Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ. 25, 486–541 (2018).29362479 10.1038/s41418-017-0012-4
30. Shen Y Cross-talk between cuproptosis and ferroptosis regulators defines the tumor microenvironment for the prediction of prognosis and therapies in lung adenocarcinoma Front. Immunol. 2022 13 1029092 10.3389/fimmu.2022.1029092 36733399
Shen, Y. et al. Cross-talk between cuproptosis and ferroptosis regulators defines the tumor microenvironment for the prediction of prognosis and therapies in lung adenocarcinoma. Front. Immunol. 13, 1029092 (2022).36733399 10.3389/fimmu.2022.1029092
31. Liu Z Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer Nat. Commun. 2022 13 816 10.1038/s41467-022-28421-6 35145098
Liu, Z. et al. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat. Commun. 13, 816 (2022).35145098 10.1038/s41467-022-28421-6
32. Xu J High PPT1 expression predicts poor clinical outcome and PPT1 inhibitor DC661 enhances sorafenib sensitivity in hepatocellular carcinoma Cancer Cell Int. 2022 22 115 10.1186/s12935-022-02508-y 35277179
Xu, J. et al. High PPT1 expression predicts poor clinical outcome and PPT1 inhibitor DC661 enhances sorafenib sensitivity in hepatocellular carcinoma. Cancer Cell Int. 22, 115 (2022).35277179 10.1186/s12935-022-02508-y
33. Sharma G PPT1 inhibition enhances the antitumor activity of anti-PD-1 antibody in melanoma JCI Insight 2020 5 e133225 10.1172/jci.insight.133225 32780726
Sharma, G. et al. PPT1 inhibition enhances the antitumor activity of anti-PD-1 antibody in melanoma. JCI Insight 5, e133225 (2020).32780726 10.1172/jci.insight.133225
34. Li X NQO1 targeting prodrug triggers innate sensing to overcome checkpoint blockade resistance Nat. Commun. 2019 10 3251 10.1038/s41467-019-11238-1 31324798
Li, X. et al. NQO1 targeting prodrug triggers innate sensing to overcome checkpoint blockade resistance. Nat. Commun. 10, 3251 (2019).31324798 10.1038/s41467-019-11238-1
35. Silvers MA The NQO1 bioactivatable drug, β-lapachone, alters the redox state of NQO1+ pancreatic cancer cells, causing perturbation in central carbon metabolism J. Biol. Chem. 2017 292 18203 18216 10.1074/jbc.M117.813923 28916726
Silvers, M. A. et al. The NQO1 bioactivatable drug, β-lapachone, alters the redox state of NQO1+ pancreatic cancer cells, causing perturbation in central carbon metabolism. J. Biol. Chem. 292, 18203–18216 (2017).28916726 10.1074/jbc.M117.813923
36. Di Biase S Fasting-mimicking diet reduces HO-1 to promote T cell-mediated tumor cytotoxicity Cancer Cell 2016 30 136 146 10.1016/j.ccell.2016.06.005 27411588
Di Biase, S. et al. Fasting-mimicking diet reduces HO-1 to promote T cell-mediated tumor cytotoxicity. Cancer Cell 30, 136–146 (2016).27411588 10.1016/j.ccell.2016.06.005
37. Tan P MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis as a novel prognostic biomarker associated with immune cell infiltration in KIRC Mol. Ther. Oncolytics 2021 22 307 325 10.1016/j.omto.2021.07.008 34553021
Tan, P. et al. MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis as a novel prognostic biomarker associated with immune cell infiltration in KIRC. Mol. Ther. Oncolytics 22, 307–325 (2021).34553021 10.1016/j.omto.2021.07.008
38. Altevogt P Doberstein K Fogel M L1CAM in human cancer Int. J. Cancer 2016 138 1565 1576 10.1002/ijc.29658 26111503
Altevogt, P., Doberstein, K. & Fogel, M. L1CAM in human cancer. Int. J. Cancer 138, 1565–1576 (2016).26111503 10.1002/ijc.29658
39. Amin MB The Eighth Edition AJCC Cancer Staging Manual: continuing to build a bridge from a population-based to a more “personalized” approach to cancer staging CA Cancer J. Clin. 2017 67 93 99 10.3322/caac.21388 28094848
Amin, M. B. et al. The Eighth Edition AJCC Cancer Staging Manual: continuing to build a bridge from a population-based to a more “personalized” approach to cancer staging. CA Cancer J. Clin. 67, 93–99 (2017).28094848 10.3322/caac.21388
40. Sammut SJ Multi-omic machine learning predictor of breast cancer therapy response Nature 2022 601 623 629 10.1038/s41586-021-04278-5 34875674
Sammut, S. J. et al. Multi-omic machine learning predictor of breast cancer therapy response. Nature 601, 623–629 (2022).34875674 10.1038/s41586-021-04278-5
41. Hänggi K Ruffell B Cell death, therapeutics, and the immune response in cancer Trends Cancer 2023 9 381 396 10.1016/j.trecan.2023.02.001 36841748
Hänggi, K. & Ruffell, B. Cell death, therapeutics, and the immune response in cancer. Trends Cancer 9, 381–396 (2023).36841748 10.1016/j.trecan.2023.02.001
42. Tan J TREM2(+) macrophages suppress CD8(+) T-cell infiltration after transarterial chemoembolisation in hepatocellular carcinoma J. Hepatol. 2023 79 126 140 10.1016/j.jhep.2023.02.032 36889359
Tan, J. et al. TREM2(+) macrophages suppress CD8(+) T-cell infiltration after transarterial chemoembolisation in hepatocellular carcinoma. J. Hepatol. 79, 126–140 (2023).36889359 10.1016/j.jhep.2023.02.032
43. Hong M YM155 inhibits topoisomerase function Anticancer Drugs 2017 28 142 152 10.1097/CAD.0000000000000441 27754993
Hong, M. et al. YM155 inhibits topoisomerase function. Anticancer Drugs 28, 142–152 (2017).27754993 10.1097/CAD.0000000000000441
44. Zabludoff SD AZD7762, a novel checkpoint kinase inhibitor, drives checkpoint abrogation and potentiates DNA-targeted therapies Mol. Cancer Ther. 2008 7 2955 2966 10.1158/1535-7163.MCT-08-0492 18790776
Zabludoff, S. D. et al. AZD7762, a novel checkpoint kinase inhibitor, drives checkpoint abrogation and potentiates DNA-targeted therapies. Mol. Cancer Ther. 7, 2955–2966 (2008).18790776 10.1158/1535-7163.MCT-08-0492
45. Liberzon A Molecular signatures database (MSigDB) 3.0 Bioinformatics 2011 27 1739 1740 10.1093/bioinformatics/btr260 21546393
Liberzon, A. et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 27, 1739–1740 (2011).21546393 10.1093/bioinformatics/btr260
46. Zou Y Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery Int. J. Surg. 2022 107 106936 10.1016/j.ijsu.2022.106936 36341760
Zou, Y. et al. Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int. J. Surg. 107, 106936 (2022).36341760 10.1016/j.ijsu.2022.106936
47. Hänzelmann S Castelo R Guinney J GSVA: gene set variation analysis for microarray and RNA-seq data BMC Bioinform. 2013 14 7 10.1186/1471-2105-14-7
Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 14, 7 (2013).10.1186/1471-2105-14-7
48. Becht E Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression Genome Biol. 2016 17 218 10.1186/s13059-016-1070-5 27765066
Becht, E. et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17, 218 (2016).27765066 10.1186/s13059-016-1070-5
49. Aran D Hu Z Butte AJ xCell: digitally portraying the tissue cellular heterogeneity landscape Genome Biol. 2017 18 220 10.1186/s13059-017-1349-1 29141660
Aran, D., Hu, Z. & Butte, A. J. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18, 220 (2017).29141660 10.1186/s13059-017-1349-1
50. Li B Comprehensive analyses of tumor immunity: implications for cancer immunotherapy Genome Biol. 2016 17 174 10.1186/s13059-016-1028-7 27549193
Li, B. et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 17, 174 (2016).27549193 10.1186/s13059-016-1028-7
51. Newman AM Robust enumeration of cell subsets from tissue expression profiles Nat. Methods 2015 12 453 457 10.1038/nmeth.3337 25822800
Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457 (2015).25822800 10.1038/nmeth.3337
52. Finotello F Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data Genome Med. 2019 11 34 10.1186/s13073-019-0638-6 31126321
Finotello, F. et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med. 11, 34 (2019).31126321 10.1186/s13073-019-0638-6
53. Racle J Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data Elife 2017 6 e26476 10.7554/eLife.26476 29130882
Racle, J. et al. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife 6, e26476 (2017).29130882 10.7554/eLife.26476
54. Wilkerson MD Hayes DN ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking Bioinformatics 2010 26 1572 1573 10.1093/bioinformatics/btq170 20427518
Wilkerson, M. D. & Hayes, D. N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573 (2010).20427518 10.1093/bioinformatics/btq170
55. Șenbabaoğlu Y Michailidis G Li JZ Critical limitations of consensus clustering in class discovery Sci. Rep. 2014 4 6207 10.1038/srep06207 25158761
Șenbabaoğlu, Y., Michailidis, G. & Li, J. Z. Critical limitations of consensus clustering in class discovery. Sci. Rep. 4, 6207 (2014).25158761 10.1038/srep06207
56. John CR M3C: Monte Carlo reference-based consensus clustering Sci. Rep. 2020 10 1816 10.1038/s41598-020-58766-1 32020004
John, C. R. et al. M3C: Monte Carlo reference-based consensus clustering. Sci. Rep. 10, 1816 (2020).32020004 10.1038/s41598-020-58766-1
57. Mayakonda A Maftools: efficient and comprehensive analysis of somatic variants in cancer Genome Res. 2018 28 1747 1756 10.1101/gr.239244.118 30341162
Mayakonda, A. et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28, 1747–1756 (2018).30341162 10.1101/gr.239244.118
58. Wu T clusterProfiler 4.0: a universal enrichment tool for interpreting omics data Innovation 2021 2 100141 34557778
Wu, T. et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation 2, 100141 (2021).34557778
59. Zhou Y Metascape provides a biologist-oriented resource for the analysis of systems-level datasets Nat. Commun. 2019 10 1523 10.1038/s41467-019-09234-6 30944313
Zhou, Y. et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10, 1523 (2019).30944313 10.1038/s41467-019-09234-6
60. Maeser D Gruener RF Huang RS oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data Brief. Bioinform. 2021 22 bbab260 10.1093/bib/bbab260 34260682
Maeser, D., Gruener, R. F. & Huang, R. S. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief. Bioinform. 22, bbab260 (2021).34260682 10.1093/bib/bbab260
